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1 .  Research  Summary 


This  research  tackles  the  problems  of  data  association  and  state 
estimation  in  the  multitarget  tracking  of  non-maneuvering  as  well  as 
maneuvering  targets  in  clutter.  The  major  contributions  are  the  development 
of  a  depth-first  search  (DFS)  algorithm  for  data  association,  three  algorithms 
which  are  suitable  for  implementation  in  multiprocessor  systems  for  fast 
computation  of  the  a  posteriori  probabilities  of  the  origins  of  measurements 
in  the  joint  probabilistic  data  association  filter  (JPDAF) ,  and  a  joint 
probabilistic  data  association  fixed-lag  smoother  (JPDAS)  for  improving  the 
accuracy  of  state  estimation. 

In  the  DFS  algorithm,  the  problem  of  data  association  is  modelled  as  an 
exhaustive  search  problem.  Based  on  this  model,  a  specialized  DFS  approach  is 
proposed  for  efficiently  generating  the  data  association  hypotheses  and 
computing  rapidly  the  conditional  probabilities  of  the  hypotheses  for  data 
association  in  the  JPDAF.  In  the  three  algorithms  developed  for 
implementation  on  multiprocessor  systems,  the  a  posteriori  probability  of  the 
origin  of  each  measurement  in  the  JPDAF  is  decomposed  into  two  parts.  The 
computation  of  one  part  becomes  trivial  and  the  different  ways  for  computing 
the  other  part  lead  to  the  development  of  the  three  algorithms.  The 
computational  costs  and  memory  requirements  of  the  above  algorithms  are 
analyzed  in  the  worst  case  as  well  as  in  the  average  case. 

A  comprehensive  analysis  reveals  several  drawbacks  of  the  neural  network 
probabilistic  data  association  (NPDA)  algorithm.  The  NPDA  is  an  application 
of  the  Hopfield  neural  network  to  the  data  association  problem. 

In  order  to  improve  the  accuracy  of  state  estimation  for  multitarget 
tracking  in  clutter,  an  algorithm  is  developed  by  making  use  of  the  joint 
probabilistic  data  association  fixed-lag  smoothing  techniques.  It  is  shown 
that  a  significant  improvement  in  the  accuracy  of  state  estimation  of  both 
nonmaneuvering  and  maneuvering  targets  may  be  achieved  by  introducing  a  time 
lag  of  one  or  two  sampling  periods  between  the  instants  of  estimation  and  the 
latest  measurement.  Finally,  a  computer  simulation  system  is  developed  to 
demonstrate  the  effectiveness  of  the  developed  algorithms. 


2 .  Research  Description 


This  research  is  mainly  focused  on  data  association  and  target  state 
estimation  in  multitarget  tracking  in  a  cluttered  environment.  The  objectives 
are  to  design  algorithms  for  efficiently  generating  the  data  association 
hypotheses,  computing  rapidly  the  a  posteriori  probabilities,  's,  and 
improving  the  accuracy  of  target  state  estimation.  The  major  contributions  of 
this  research  are  the  development  of 

•  the  depth-first  search  (DFS)  algorithm  for  data  association, 

•  three  algorithms  for  computing  the  a  posteriori  probabilities  using 
multiprocessor  systems,  and 

•  the  JPDA  fixed-lag  smoothing  (JPDAS)  algorithm. 

In  the  DFS  algorithm,  the  problem  of  data  association  is  identified  as 
an  exhaustive  search  problem  in  general.  Subsequently,  a  mathematical  model 
is  proposed  for  the  problem  of  data  association  in  the  JPDAF .  Based  on  the 
model,  a  specialized  DFS  approach  is  developed  for  the  fast  generation  of  data 
association  hypotheses  and  for  the  computation  of  the  conditional 
probabilities  of  the  hypotheses  in  the  JPDAF.  The  computational  complexity  of 
the  algorithm  is  analyzed  in  terms  of  the  number  of  multiplications  and 
additions  required  in  the  computation  of  the  a  posteriori  probabilities, 
pj's,  in  the  worst  case  as  well  as  in  the  average  case.  Although  the  DFS 
algorithm  requires  more  multiplications  and  additions  that  the  fast  JPDA 
algorithm  [1]  in  the  worst  case,  it  is  much  more  efficient  than  the  fast  JPDA 
algorithm  in  the  average  case.  In  addition,  the  DFS  algorithm  requires  less 
memory  than  the  fast  JPDA  algorithm.  Another  advantage  of  the  DFS  algorithm 
is  that  it  could  be  easily  extended  to  generate  the  data  association 
hypotheses  in  the  measurement-oriented  approach  [2] . 

Three  algorithms  which  are  suitable  for  implementation  in  a 
multiprocessor  system  are  also  developed.  In  the  three  algorithms,  the 
computation  of  the  a  posteriori  probabilities,  pr  s  is  not  based  on  the 
generation  of  the  data  association  hypotheses  like  in  the  DFS  algorithm.  Each  pj 
in  the  three  algorithms  is  decomposed  into  two  parts.  The  computation  of  one 
part  is  trivial  and  the  various  ways  of  computing  the  other  part  (due  to 
interference  from  other  targets)  lead  to  the  development  of  the  three 
algorithms . 

In  the  first  algorithm,  the  computation  of  the  interference  part  of  pj 
is  performed  recursively  in  a  top-to-bottom  mode.  In  the  worst  case,  this 
recursive  algorithm  requires  more  multiplications  and  additions  than  the  fast 
JPDA  algorithm  [1].  However,  in  the  average  case,  the  recursive  algorithm  is 
expected  to  perform  much  better  than  the  fast  JPDA  algorithm.  Furthermore, 
the  recursive  algorithm  requires  less  storage  space  than  the  fast  JPDA 
algorithm.  In  comparison  with  the  DFS  algorithm  developed  earlier  in  this 
research,  the  recursive  algorithm  requires  less  multiplications  and  additions 
but  more  storage  space  than  the  DFS  algorithm.  The  most  important  feature  of 
the  recursive  algorithm  is  that  it  can  be  implemented  on  a  multiprocessor 
system,  which  could  speed  up  the  computation  of  the  pi'3  significantly.  On 
the  other  hand,  the  recursive  algorithm  is  not  as  suitable  as  the  DFS 
algorithm  for  extension  to  the  measurement-oriented  approach  since  it  is  not 
based  on  the  generation  of  data  association  hypotheses. 

The  second  algorithm  is  a  nonrecursive  algorithm.  Since  the 
interference  part  of  pj  is  computed  by  polynomial  multiplication,  the 
computational  cost  is  drastically  reduced  with  only  a  slight  increase  in 
memory  when  the  density  of  targets  is  not  very  high.  The  nonrecursive 
algorithm  is  especially  suitable  for  a  small  tracking  system  where  the 
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computational  capability  is  limited. 


In  the  last  algorithm,  only  the  interference  from  the  neighboring 
targets  of  target  t  is  considered  in  the  computation  of  the  (Jj's.  The 
computational  cost  and  memory  requirement  are  only  marginally  increased  in 
this  approximating  algorithm  when  it  is  compared  with  the  PDAF  [3] .  It  is 
apparent  from  the  computer  simulation  that  performance  of  the  approximating 
algorithm  is  quite  close  to  that  of  the  original  JPDAF  in  most  of  the 
scenarios  presented  in  the  simulation.  Generally  speaking,  the  approximating 
algorithm  is  proposed  for  applications  where  both  the  computational  efficiency 
and  the  memory  requirement  are  critical. 

Another  issue,  which  is  quite  different  from  data  association  as 
discussed  above,  in  multitarget  tracking  is  brought  to  attention.  In  the 
original  JPDAF,  the  accuracy  of  the  target  state  estimation  is  affected  not 
only  by  the  measurement  noise  but  also  by  the  uncertainty  in  the  origins  of 
the  measurements.  To  improve  the  accuracy  of  the  target  state  estimation,  a 
JPDAS  (JFDA  smoothing)  algorithm  is  proposed.  Some  computer  simulation 
results  have  been  presented  to  illustrate  the  improvement  in  tracking  accuracy 
achieved  by  introducing  a  small  time  delay  between  the  instants  of  estimation 
and  the  latest  measurement.  The  experimental  results  show  that  the  JPDAS 
algorithm  works  well  for  tracking  multiple  targets  in  a  cluttered  environment. 
However,  the  introduction  of  the  time  lag  causes  an  increase  in  the 
computational  and  storage  requirements.  Fortunately,  a  small  lag  produces  a 
significant  improvement  in  the  estimation  accuracy  so  that  the  increase  in  the 
computational  burden  is  kept  to  a  minimum. 

In  addition  to  the  development  of  the  algorithms  for  data  association 
and  target  state  estimation,  a  comprehensive  analysis  of  the  neural  network 
solution  to  the  data  association  problem  in  [4]  and  [5]  is  also  provided.  The 
analysis  reveals  that  the  Hopfield  neural  networks  developed  in  [4]  and  [5] 
have  improper  energy  functions.  This  resulted  from  misinterpretations  of  the 
properties  of  the  JPDAF  which  the  networks  were  designed  to  emulate  and  the 
improper  selections  of  the  constant  coefficients  in  the  energy  functions.  The 
computer  simulation  system  developed  for  multitarget  tracking  in  a  cluttered 
environment  is  available  along  with  the  simulation  results  from  the  system. 
This  simulation  system  consists  of  six  modules,  spanning  system  specification, 
measurement  generation,  clustering,  data  association,  target  state  estimation, 
and  run-time  monitoring. 
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3 .  The  Most  Significant  Research  Accomplishment 

The  most  significant  research  accomplishment  is  the  incorporation  of 
image-based  tracking  with  errors  in  measurement  as  well  as  errors  in 
registration,  which  can  be  coupled  with  the  efficient  data  association 
algorithm  developed  for  the  target-oriented  approach  (and  which  is  suitable 
for  generalization  to  the  measurement-oriented,  trak-oriented,  and  multiple 
correlation  approaches)  for  fast  and  accurate  tracking  of  nonmaneuvering  as 
well  as  maneuvering  multitargets  in  clutter.  The  image  processing  algorithms 
are  necessary  to  identify  targets  in  the  presence  of  anticipated  background 
noise  (including  earth,  lunar,  star  backgrounds,  complicated  spacecraft 
structures,  specular  reflections)  and  blur  (including  atmospheric  turbulence, 
motion,  optical  aberration) .  The  areas  of  application  include  not  only 
surveillance,  but  also  autonomous  rendezvous  and  capture,  autonomous  planetary 
landing  systems,  and  satellite  servicing.  The  sensors  could  be  passive  (video 
cameras,  stereo  vision  (3-D  sensor)),  as  well  as  active  (laser  radar)  and  the 
benefits  of  image-based  tracking  systems  could  include  economy  (cheaper 
tracking  systems)  in  addition  to  speed  and  accuracy. 


-4- 


4 .  Research  Personnel 


(a)  Dr.  N.  K.  Bose,  HRB-Systems  Professor  of  Electrical  Engineering  and 
Director  of  The  Spatial  and  Temporal  Signal  Processing  Center  at  The 
Pennsylvania  State  University,  University  Park,  PA  16802. 

Was  appointed  Principal  Investigator  of  the  project,  originally 
under  the  supervision  of  Professor  A.  K.  Mahalanabis,  who  passed  away  on 
July  31,  1988.  Dr.  Bose  supervised  the  original  project,  whose  funding 
expired  in  December  1988,  and  prepared  a  new  proposal  entitled, 
"Multitarget  Tracking  in  Clutter:  New  Approaches,"  which  was  approved  for 
funding  for  a  3-year  period,  effective  March  1989.  A  no-cost  extension 
permitted  extension  of  the  contract  (Number  N00014-86-K-0542)  to  May  15, 
1992.  During  the  period  of  the  contract.  Dr.  Bose  supervised  the 
research,  provided  new  ideas  for  implementing  data  association 
efficiently,  and  prepared  manuscripts  for  submissionto  and  subsequent 
publication  in  reviewed  journals.  Dr.  Bose  also  interacted  with  Professor 
Leon  H.  Sibul,  of  the  Applied  Research  Laboratory  at  The  Pennsylvania 
State  University,  and  other  scientists  working  on  multitarget  tracking, 
especially  in  the  SDIO/IST  Program.  He  regularly  reported  the  progress  of 
research  at  the  annual  SDIO/IST  Workshop  at  Arlington,  Virginia,  where 
other  principal  investigators  also  came  to  present  their  ideas  and  share 
research  results. 

(b)  Bin  Zhou,  post-doctoral  scholar  at  The  Center  for  Multivariate  Analysis 
and  The  Spatial  and  Temporal  Signal  Processing  Center  at  The  Pennsylvania 
State  University,  University  Park,  PA  16802. 

Was  continuously  supported  as  a  graduate  research  assistant  until 
the  completion  of  his  Ph.D.  dissertation  in  late  1991.  This  dissertation 
is  on  very  high  demand  and  requests  for  copies  of  the  dissertation  have 
been  received  from  scientists  from  all  over  the  world,  including 
Indonesia,  Italy,  and  Sweden.  Bin  Zhou  has  developed  a  complete  computer 
simulation  system  to  simulate  multiple  target  tracking  in  a  cluttered 
environment . 
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5.  Publications 


(From  previous  ONR  Contract  No.  N0014-86— K-0542-P004,  July  1,  1989  to  May  15, 
1992) 


1.  N.  K.  Bose,  H.  C.  Kim,  and  H.  M.  Valenzuela,  "Recursive  implementation  of 
total  least  squares  algorithm  for  image  reconstruction,"  IEEE  Transactions 
on  Image  Processing,  under  review. 

2.  N.  K.  Bose  and  L.  H.  Sibul,  "Multidimensional  Signal  Processing:  Sensor 
Array  Processing,  "  invited  chapter  in  Handbook  of  Electrical  Engineering, 
ed.  R.  C.  Dorf,  CRC  Press,  scheduled  to  appear  in  1993. 

3.  B.  Zhou  and  N.  K.  Bose,  "A  Comprehensive  Analysis  of  the  'Neural  Solution 
to  the  Multitarget  Tracking  Data  Association  Problem',"  IEEE  Trans,  on 
Aerospace  and  Electronic  Systems,  accepted  for  publication  and  scheduled 
to  appear  in  late  1992. 

4.  B.  Zhou  and  N.  K.  Bose,  "Multitarget  Tracking  in  Clutter:  Fast  Algorithms 
for  Data  Association,"  IEEE  Trans,  on  Aerospace  and  Electronic  Systems, 
accepted  for  publication  and  scheduled  to  appear  in  mid-1993. 

5.  A.  K.  Mahalanabis,  B.  Zhou,  and  N.  K.  Bose,  "Improved  Multi-Target 
Tracking  in  Clutter  by  PDA  Smoothing,"  IEEE  Trans,  on  Aerospace  and 
Electronic  Systems,  Vol.  26,  January  1990,  pp.  113-121. 

6.  A.  K.  Mahalanabis  and  B.  Zhou,  "A  Joint  Probabilistic  Data  Association 
Smoothing  Algorithm  for  Multitarget  Tracking,"  the  1988  American  Control 
Conference,  Atlanta,  GA,  June  1988,  pp.  360-365. 


Ph.D.  Dissertation 

B.  Zhou,  "Multitarget  Tracking  in  a  Cluttered  Environment,"  Ph.D. 
Dissertation,  Department  of  Electrical  and  Computer  Engineering,  The 
Pennsylvania  State  University,  University  Park,  PA,  1991. 
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6 .  Selected  preprints 
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Image  Processing 

Recursive  Implementation  of  Total  Least  Squares  Algorithm 
for  Image  Reconstruction 
From  Noisy,  Undersampled  Multiframes  * 

N.  K.  Bose,  H.  C.  Kim,  and  H.  M.  Valenzuela 

Department  of  Electrical  Engineering 
The  Spatial  and  Temporal  Signal  Processing  Center 
The  Pennsylvania  State  University 
University  Park,  PA  16802 
Tel:(814)  865-3912 
Fax:(814)  865-7065 

ABSTRACT t 

It  is  shown  how  the  efficient  recursive  total  least  squares  algorithm  recently 
developed  by  C.  E.  Davila  for  real  data  can  be  applied  to  image  reconstruction  from 
noisy,  undersampled  multiframes  when  the  displacement  of  each  frame  relative  to  a 
reference  frame  is  not  accurately  known.  To  do  this,  the  complex- valued  image  data 
in  the  wavenumber  domain  is  transformed  into  an  equivalent  real  data  problem  to 
which  Davila’s  algorithm  is  successfully  applied.  Two  detailed  illustrative  examples 
are  provided  in  support  of  the  procedure. 
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1.  Introduction 


In  many  applications,  multiple  low  resolution  images  of  an  object  are  available, 
but  a  high  resolution  image  is  desired.  For  example,  a  camera  aboard  a  satellite  is 
used  to  take  several  pictures  of  an  area  on  the  ground  and  to  transmit  a  sequence  of 
undersampled  low  resolution  frames,  which  are  shifted  from  each  other.  Futhermore, 
those  frames  are  often  degraded  versions  of  the  original  scene  due  to  blur  and  noise. 
Through  the  task  of  image  registration,  the  displacement  of  each  frame  relative 
to  an  arbitrarily  chosen  reference  frame  is  measured.  The  type  of  blur  and  the 
characteristic  of  noise  can  be  determined  from  the  understanding  of  the  various 
physical  processes  involved  in  the  image  formation.  Subsequently,  interpolation, 
deblurring,  and  filtering  are  required  to  reconstruct  a  high  resolution  image. 

S.  P.  Kim,  N.  K.  Bose,  and  H.  M.  Valenzuela  recently  developed  an  efficient 
algorithm  to  reconstruct  a  high  resolution  noise-free  image  from  undersampled 
low  resolution  noisy  multiframes  [1].  This  algorithm  implements,  simultaneously, 
the  tasks  of  interpolation  and  noise  filtering  of  the  input  images  by  applying 
the  recursive  least  squares  (RLS)  sequential  estimation  theory  in  the  wavenumber 
domain. 

In  [1],  the  displacements  of  the  images  are  assumed  to  be  known.  Actually, 
errors  occur  in  the  estimation  of  the  relative  displacement  during  image  registration. 
To  combat  this  problem,  we  apply  the  Total  Least  Squares  (TLS)  method.  This 
method  is  known  to  be  very  useful  for  improving  the  solution  accuracy  when  errors 
are  present  not  only  in  the  observation  but  also  in  the  measurement  matrix  [2].  A 
Recursive  Total  Least  Squares  (RTLS)  algorithm  was  developed  by  C.  E.  Davila  [3] 
in  the  context  of  the  adaptive  filtering  problem. 
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In  Section  2,  the  image  interpolation  model  which  is  formulated  by  using  the 
aliasing  property  of  the  Discrete  Fourier  Transform  (DFT)  is  briefly  described. 
The  algorithm  for  high  resolution  image  reconstruction,  described  in  Section  3,  is 
obtained  by  transforming  the  original  complex-data  problem  to  an  equivalent  real 
data  problem  to  which  the  RTLS  algorithm  in  [3]  is  applied.  In  Section  4,  the 
performance  of  the  algorithm  is  verified  by  computer  simulations. 

2.  Reconstruction  Model  of  High  Resolution  Images 

The  A:-th  (M  x  N )  undersampled  observed  image  frame,  fk{hj),  *  =  0,1,..., 
M  —  1,  j  =  0, 1, . . . ,  N  —  1,  of  the  noise-free  original  continuous  image,  f(x,  y),  can 
be  written  as 

fk(i,j)  =  f(iTx  +  6xk,jTy  +  6yk),  fc  =  l,2,...,p  (1) 

where  6xk  and  Sy k  are  the  estimated  shifts  and  Tx  and  Ty  are  the  sampling  periods 
along  the  x  and  y  axes,  respectively  [1].  We  assume  that  the  Fourier  transform, 
Fc(u,  v )  of  the  original  continuous  image  f(x,y)  is  approximated  by  the  bandlimited 
constraint 

|Fc(u,v)|  =  0,  |u|  >  Lxwz  and  |v|  >  Lywy  (2) 

for  some  finite  integers  Lx  and  Ly  where  wx  =  (27r/rr)  and  wy  =  (2ir/Ty). 

The  multiframe  image  restoration  model  given  in  [1]  can  be  written  in  the  form 

Zk  =  Y{F  +  Nk  (3) 

where  the  superscript  t  denotes  the  transpose  operation,  Zk  is  the  value  at  the 
wavenumber  point  (m,  n)  of  the  DFT  of  the  fc-th  noisy  undersampled  frame,  iV*  is 
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the  corresponding  noise  term  in  the  wavenumber  domain,  F_  is  a  vector  containing 
the  p  =  4 LxLy  interpolated  components  at  wavenumber  point  (m,n)  denoted  by 

F=  [Fmn(l),  Fmn(2),...,Fmn(p)} 

a  , »  «  -  *  (4) 

t/i  i  fi  >  •  •  •  ?  fp] 

and  the  coefficient  vector  Yk  is  defined  as 

Yk  =  [^fc,u  <t>k, 2,  ■  •  •  ,<f>k,p]  (5) 

where 

fc.- = fkexr{j2,r{SMitk + k)+l’k(ik + ^)}}  (6) 

after  defining  the  lexicographical  ordering  associated  with  row-wise  scanning  by 

lx  —  {r  -  l)mod{2Lx)  -  Lx  and  ly  =  [(r  -  1)/(2LX)\  -  Ly.  (7) 

If  there  are  J  frames  available(J  >  4 LxLy),  then  eqn.  (3)  may  be  used  to  obtain 
the  matrix  equation 

Zj  =  *jF  +  Nj  (8) 

where 

±j  =  [YuY2,...,Yj}t 
Zj  =  [Zi,z2,...  ,ZjY 

By  solving  eqn.  (8),  we  obtain  the  high  resolution  samples  at  the  generic 
wavenumber  point  (m,n).  The  interpolating  factors  are,  respectively,  2 Lx  and  2 Ly 
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along  the  x  and  y  axes.  The  inverse  DFT  oi  the  interpolated  image  gives  the  high 
resolution  image. 


3.  Reconstruction  of  High  ResolutionTmage  by  RTLS  Algorithm 


The  TLS  algorithm  provides  a  solution  when  there  are  errors  in  both  the 
observation  vector  and  in  the  measurement  matrix.  Prior  to  [3],  the  solution  was  not 
generated  recursively.  In  [3],  an  algorithm  was  developed  for  the  linear  regression 
problem  to  produce  unbiased  filter  coefficients  for  FIR  adaptive  filters.  To  be  able 
to  apply  the  result  in  [3]  to  the  image  reconstruction  problem  under  consideration 
here  (with  both  observation  noise  and  inter-frame  displacement  error  present)  it  is 
necessary  to  follow  the  procedure  described  below. 

To  account  for  the  errors  in  estimation  of  the  shifts  during  the  registration 
phase,  we  define, 


&zi t  —  6zk  +  &6xk 

6yk  —  6yk  d"  A5yjfc, 


(9) 


where  8xk  and  8yk  are  the  actual  shifts  and  A 6xk,  A6yk  are  the  respective 
displacement  errors  along  the  x  and  y  axes  in  the  registration  process.  Substituting 
eqn.  (9)  into  eqn.  (6),  we  obtain 


4>kr  =  ^T^exp{j27r{(6xfc-b A6r*)(-^r  + A-)-t-(5yAr  +  A^jfe)(^r  +  ^-)}|.  (10) 

For  small  magnitudes  of  A 6xk  and  A 6yk,  <j>kr  in  eqn.  (10)  may  be  approximated  by 


<f>kr  —  <f>kr  +  ekr 


(11) 
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where  <pkr,  defined  in  eqn.  (6),  is  expressible  as 


<f>kr  =  <f>kr 


A<5rJfe  =  A6y  fc  =0 


and 


for 


Ckr  -  j2ir{A6xk{j^r  +  ~)  +  A +  ^L)}d'fcr 


NTy  Ty 


2n{ASxk{  MT *  +  +  AM  NT^  +  ^  )} 


«  1. 


(12) 


(13) 


Since  the  maximum  values  of  m  and  n  are,  respectively,  M  —  1  and  N  —  1,  the 
inequality  in  eqn.  (13)  has  only  to  be  checked  for  the  case  when  m  =  M  —  1  and 
n  =  N  —  1.  Substituting  eqn.  (7)  in  eqn.  (13),  simple  sufficient  conditions  for  the 
bound  in  eqn.  (13)  to  hold  are 


and 


\A6zk\  « 


|A^i  < 


IB. 


4t r(Lx  +  1)  2  (Lx  +  1) 


(14a) 


IB, 


47r(Zy  +1)  2(  Ij  y  +  1) 

Using  eqns.  (8)  and  (11)  and  after  replacing  the  subscript  J  by  k ,  we  get 


(146) 


Zk  =  [*k  +  Ek}F  +  Nk 


(15) 


where  the  ( i,j)th  element  of  the  matrix  Ek  is  of  eqn.  (12).  The  matrix  Ek  in 
eqn.  (15)  is  basically  a  perturbation  of  the  coefficient  matrix  produced  by  the 
delay  estimation  errors.  We  now  consider  the  TLS  problem  formulated  as 


minimize  \\{Nk :  £*]|| 
subject  to  -  N_k  =  -f  Ek]F 


(16) 
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where  ||  •  ||  denotes  the  Frobenius  norm 

*  J 

on  the  matrix  B  =  [6,-j].  Define 

A  • 

it  =  =  4*1 

(17) 

fZk  =  HSU  '■  £»] 

(18) 

(19) 


(20) 


and 

«  =  (Jp)  =[1,-/1, -A, 

The  TLS  problem  in  eqn.  (16)  can  be  restated  as: 

minimize  |(|Ffc|| 
subject  to  +  Kk)qk  =  0. 

From  what  is  known  in  [4] ,  the  above  constrained  minimization  problem  can  be 
associated  with  the  equivalent  minimization  problem 

A  yj 

n(qk)  =  minlljFji.il2  =  min  qkMt—k9k  =  min  (21) 

9k  Qk  9k  9k 

where  the  superscript  H  denotes  the  complex  conjugate  transpose  operation, 
referred  to  as  Hermitian  conjugation.  The  Hermitian  matrix  Rk  in  eqn.  (21)  is 
defined  as 

Rk  =  if  S*  =  £  r)r)  (22) 

}=  1 

where  the  star  superscript  denotes  the  complex  conjugate  operation  and 


rk  —  [ Zk ,  <t>k  1  +£kl,  <t>k2  +  Ck2,-  •  •  ,<f>kj>  +  ffcp]*- 


(23) 
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The  minimization  problem  in  eqn.  (21)  corresponds  to  the  obtaining  of  the 
eigenvector  qk  associated  with  the  smallest  eigenvalue  of  Rk •  Given  the  previous 
eigenvector  qk- 1,  we  update  it  to  obtain  qk  from 

qk  =  qk-i  +  <xki>k,  (24) 

as  in  [5],  where  V’fc  is  a  correction  vector,  chosen  in  [3]  to  be  the  Kalman  gain  vector, 

rf>k  =  Kk  =  R;'rk  (25) 

and  the  scalar  a*  is  derived  from  the  minimization  of 

*7(9*)  =  v(Qk- 1  +  akij>k)  =  ~  -•  (26) 

In  order  to  update  the  Kalman  gain,  in  our  formulation  we  use  the  matrix  inversion 

lemma  in  the  R.LS  algorithm.  To  get  qk  recursively,  as  in  eqns.  (24)-(26),  we 

provide  a  mechanism  for  transforming  our  problem  involving  complex  variables  to 
an  equivalent  problem  involving  real  variables  where  C.  E.  Davila’s  algorithm  [3] 
directly  applies.  To  do  this,  let  3?(x)  and  3(x)  denote,  respectively,  the  real  and 
imaginary  parts  of  x.  Then  express  Rk  and  qk  as 

Rk  =  Ak+jBk  (27) 


where 

Ak  =  »(«*)  and  Bk=S(Rk), 


and 

qk  =  ck  +  jdk 


(28) 
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where 


Ck  =  &(?*)  and  dk  -  S(gjfe)- 


Gk 


“( 


Since  Rk  is  Hermitian,  therefore  Ak  =  Ak  and  B\  =  —  Bk.  Define 

Ak  Bk\ 

Ak) 

and 

Then  it  can  be  verified  that  q(qk)  in  eqn.  (26)  is  expressible  as 


vM 


qjfRkqk  _  h}kGkhk 


q*fqk  h^hk 

If  we  express  the  “Kalman  gain”  Kk  in  our  problem  as 


Kk  =  Vk  +jWk 


where 

Vk=X{Kk)  and  Wk=S{Kk), 
then  an  associated  vector  Mk  is  defined  below. 


Given  hk~\,  the  updated  vector  hk  in  eqn.  (30)  is 

hk  -  hjt-i  +  &kMk , 


where  the  scalar  0k  is  obtained  by  minimizing 

’l(hk)  =  v(hk-i  +0kMk)  = 


(29) 

(30) 

(31) 

(32) 

(33) 

(34) 

(35) 
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To  update  Gk  in  eqn.  (29),  we  need  to  update  Ak  and  Bk  from  rk.  We  define 


sk  =  9?(rjt)  and  ek  =  S(rfc).  (36) 

It  can  be  easily  shown  that 

Ak  —  Ak- 1  +  sits*  +  (37) 

Bk  =  Bk-i  A  Skel  -  ek3k-  (38) 

Substituting  eqn.  (34)  into  eqn.  (35)  and  differentiating  with  respect  to  /3,  we  obtain 
the  quadratic  equation 

a/32  +  b/3  +  c  =  0  (39) 

where 

a  =  htk_lGk<t>k<t>k<t>k  ~  <t>tkGk<f>khtk^1<t>k 
b  =  htk_1Gkhk-i4>tk<f>k  -  <t>kGk<t>k  (40) 

c  =  hk-\Gkhk-\htk_l<f>k-\  -  htk_1Gk4>k , 

as  in  [3].  It  has  been  shown  in  [5]  that  eqn.  (39)  has  only  recil  roots.  Choosing  the 
smallest  root  of  this  equation,  we  obtain  the  updated  vector  hk  from  eqn.  (34).  We 
now  summarize  the  algorithm  below. 
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MODIFIED  ALGORITHM  FOR  COMPLEX  DATA 

<t>k  —  Mk 

Ak  =  +  3k^k  +  efcefc 

Bk  =  Bk- 1  +  ske\  -  ek s\ 


Gk  — 


A*  Bk 

Bi  Ak 


a  =  h^^k^k^k^k  ~  <f>tkGk<f>khtk-i<i>k 

b  =  h\_1Gkhk-i(i>tk4>k  -  4>\Gk<i>k 

c  =  —  ^1-1  Gk(j)k 

_  ~b  —  y/b2  —  Aac 


hk  =  hk- 1  +/?jb^>* 


y/h\hk 


The  final  solution  can  be  obtained  by  converting  hk  to  qk  and 

i  _  MU.  ■  _  0 
/*~1  («).’  ’■",p 

where  (gjt  )i  is  the  i-th  element  of  the  vector  qk. 


4.  Computer  Simulation 

In  this  computer  simulated  examples,  a  high  resolution  (256  x  256)  image  is 
recursively  reconstructed  from  a  set  of  low  resolution  (128  x  128)  noisy  frames  which 
are  shifted  with  respect  to  a  reference  frame.  These  shifts  or  displacements  are  not 
accurately  known.  In  order  to  reduce  the  size  of  the  arrays  to  be  processed,  the 
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input  image  was  partitioned  into  16  nonoverlapping  sections  each  of  size  (32  x  32) 
and  the  recursive  Total  Least  Squares  (RTLS)  algorithm  described  in  Section  3  was 
independently  applied  to  each  such  section.  It  is  recalled  that  in  [1]  a  similar 
approach  was  adopted  but  instead  of  RTLS  *the  recursive  least  squares  (RLS) 
algorithm  was  used.  For  each  one  of  these  16  sections,  the  interpolation  problem 
corresponds  to  the  reconstruction  of  a  (64  x  64)  image  from  a  sequence  of  shifted 
low  resolution  (32  x  32)  noisy  input  frames  when  the  inter-frame  displacements 
are  not  accurately  known.  To  generate  the  K  —  16  shifted  low  resolution  input 
frames  required  in  the  simulation  from  the  available  data,  we  use  the  DFT-based 
interpolation  technique  described  in  [1].  After  assigning  one  of  the  input  frames  to 
be  the  reference  frame,  we  label  it  as  frame  number  1  and  the  remaining  ones  are 
sequentially  labelled  from  k  =  2  to  fc  =  16.  The  relative  shifts  of  these  frames  with 
respect  to  the  reference  one  along  the  x  and  y  axes  are  denoted,  respectively,  by 
6xk  and  6yk  for  k  —  2, 3, . . . ,  16. 

To  simulate  the  errors  A6**  and  A bvk,  k  =  2, 3, . . . ,  16,  in  the  shifts  in  eqn.  (9), 
we  use  the  expression  in  eqn.  (12)  after  assigning  to  these  displacement  errors 
uniformly  distributed  random  values  over  the  interval  [—  |,j]-  Subsequently,  an 
uniformly  distributed  zero  mean  noise  was  added  to  each  input  frame.  These  noise 
sequences  were  statistically  independent  and  the  signal-to-noise  ratio  (SNR)  was 
chosen  to  be  20dB.  The  first  simulation  example  uses  the  image  of  a  girl  given  in 
Fig.  1.  Fig.  2  shows  4  (out  of  a  total  of  16)  low  resolution  (128  x  128)  noisy  input 
frames,  each  having  a  displacement  error  with  respect  to  the  chosen  reference  frame 
(naturally,  the  reference  frame  is  defined  to  have  zero  displacement  error),  and  an 
additive  noise  with  SNR  level  of  20dB.  The  RTLS  algorithm  given  in  Section  3  was 
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used  to  recursively  obtain  new  estimates  of  the  reconstructed  high  resolution  image. 
The  last  estimate  (the  sixteenth)  of  this  sequence  is  shown  in  Fig.  3.  The  visual 
quality  of  the  reconstructed  image  in  Fig.  3  is  very  close  to  the  original  image  in 
Fig.  1. 

To  illustrate  the  intermediate  steps  of  the  reconstruction,  we  provide  in  Fig.  4 
a  generic  set  of  16  (32  x  32)  frames  corresponding  to  a  section  of  the  low  resolution 
noisy  image.  From  this  set,  a  (64  x  64)-section  image  is  obtained  by  application 
of  the  RTLS  algorithm.  The  sequence  of  estimates  leading  to  the  construction  of 
the  high  resolution  filtered  image  is  shown  in  Fig.  5.  As  expected,  the  first  few 
iterations  of  the  RTLS  algorithm  provide  poor  estimates  because  the  algorithm  is 
in  a  transient  stage.  Subsequently,  the  estimates  improve  and  the  algorithm  could 
be  terminated  even  before  the  sixteenth  iteration  without  seriously  affecting  the 
visual  quality  of  the  reconstructed  image. 

The  second  simulation  example  uses  the  image  of  an  aerial  photograph  given 
in  Fig.  6.  Fig.  7  shows  4  (out  of  a  total  of  16)  low  resolution  (128  x  128)  noisy  input 
frames,  each  having  a  displacement  error  with  respect  to  the  chosen  reference  frame, 
and  an  additive  noise  with  SNR  level  of  20dB.  The  last  estimate  (the  sixteenth) 
obtained  by  application  of  the  RTLS  algorithm  is  shown  in  Fig.  8.  Fig.  9  shows  a 
generic  set  of  16  (32  x  32)  frames  corresponding  to  a  section  of  the  low  resolution 
noisy  image.  From  this  set,  a  (64  x  64)-section  image  is  obtained  by  application  of 
the  RTLS  algorithm.  The  sequence  of  estimates  leading  to  the  construction  of  the 
high  resolution  filtered  image  for  this  section  is  shown  in  Fig.  10. 
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5.  Conclusion 


It  is  shown  how  the  total  least  squares  recursive  algorithm  for  the  real  data  FIR 
adaptive  filtering  problem  may  be  applied  to  reconstruct  a  high  resolution  filtered 
image  from  undersampled,  noisy,  multiframes,  when  the  inter-frame  displacements 
are  not  accurately  known.  This  is  done  in  the  wavenumber  domain  after 
transforming  the  complex  data  problem  to  an  equivalent  real  data  problem,  to  which 
the  algorithm  developed  in  [3]  applies.  The  procedure  developed  also  applies  to  the 
case  when  the  multiframes  are  degraded  by  linear  shift-invariant  blurs,  in  a  manner 
similar  to  that  achieved  in  [1],  where  the  least  squares  sequential  estimation  theory 
was  applied.  All  the  advantages  regarding  implementation  via  massively  parallel 
computational  architecture  apply  here  as  in  [1],  because  of  the  decoupling  resulting 
from  the  type  of  processing  carried  out  in  the  wavenumber  domain  -  a  feature  also 
in  evidence  in  a  different  but  related  context  [6]. 
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Fig.  1.  Original  image  -  a  girl 


Fig.  2.  1st  -  1th  (128  i2S)  input  frames 

(  r  1  pixel  displacement  error  and  20dB  SNR) 


Fig.  3.  The  final  (256  x  256)  reconstructed  image  by  RTLS 
(after  16  frames  are  processed.) 
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Fig.  4.  16  (32  x  32)  input  frames 

(±|  pixel  displacement  error  and  20dB  SNR) 
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Fig.  5.  (64  •  64)  reconstructed  image  by  R1LS 
(each  output  is  obtained  recursively.) 
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Fig.  9.  16  (32  >  32)  input  frames 

( —  1  pixel  displacement  error  and  20dB  SNR) 


Fig.  10.  (64  64)  reconstructed  image  by  RTLS 

(each  output  is  obtained  recursively.) 


To  appear  in  The  Electrical  Engineering  Handbook,  January  1993. 


MULTIDIMENSIONAL  SIGNAL  PROCESSING:  SENSOR  ARRAY  PROCESSING 

N.  K .  Bose  and  L.  H.  Sibul,  The  Pennsylvania  State  University 

Introduction 

Multidimensional  signal  processing  tools  apply  to  aperture  and  sensor 
array  processing.  Planar  sensor  arrays  can  be  considered  to  be  sampled 
apertures.  Three  dimensional  or  volumetric  arrays  can  be  viewed  as 
multidimensional  spatial  filters.  Therefore,  the  topics  of  sensor  array 
processing,  aperture  processing,  and  multidimensional  signal  processing  can  be 
studied  under  a  unified  format.  The  basic  function  of  the  receiving  array  is 
transduction  of  propagating  waves  in  the  medium  into  electrical  signals. 
Propagating  waves  are  fundamental  in  radar,  communication,  optics,  sonar,  and 
geophysics.  In  electromagnetic  applications,  basic  transducers  are  antennas 
and  arrays  of  antennas.  A  large  body  of  literature  that  exists  on  antennas 
and  antenna  arrays  can  be  exploited  in  the  areas  of  aperture  and  sensor  array 
processing.  Much  of  the  antenna  literature  deals  with  transmitting  antennas 
and  their  radiation  patterns.  Due  to  the  reciprocity  of  transmitting  and 
receiving  transducers,  key  results  that  have  been  developed  for  transmitters 
can  be  used  for  analysis  cf  receiver  aperture  and/or  array  processing. 
Transmitting  transducers  radiate  energy  in  desired  directions,  whereas 
receiving  apertures/arrays  act  as  spatial  filters  that  emphasize  signals  from 
a  desired  look  direction  while  discriminating  against  interferences  from  other 
directions.  The  spatial  filter  wavenumber  response  is  called  the  receiver 
beampattern.  Transmitting  apertures  are  characterized  by  their  radiation 
patterns . 

Conventional  beamforming  deals  with  the  design  of  fixed  beampatterns  for 
given  specifications.  Optimum  beamforming  is  the  design  of  beampatterns  to 
meet  a  specified  optimization  criterion.  It  can  be  compared  to  optimum 
filtering,  detection,  and  estimation.  Adaptive  beamformers  sense  their 
operating  environment  (for  example,  noise  covariance  matrix)  and  adjust 
beamformer  parameters  so  that  their  performance  is  optimized  [Monzingo  and 
Miller,  1980].  Adaptive  beamformers  can  be  compared  with  adaptive  filters. 

Multidimensional  signal  processing  techniques  have  found  wide 
application  in  seismology,  where  a  group  of  identical  seismometers,  called 
seismic  arrays,  are  used  for  event  location,  studies  of  the  earth's 
sedimentation  structure,  and  separation  of  coherent  signals  from  noise,  which, 
sometimes,  may  also  propagate  coherently  across  the  array  but  with  different 
horizontal  velocities,  by  emplov^ng  velocity  filtering  [Claerbout,  1976], 
Velocity  filtering  is  performed  by  multidimensional  filters,  and  allows  also 
for  the  enhancement  of  signals  which  may  occupy  the  same  wavenumber  range  as 
noise  or  undesired  signals  do.  In  a  broader  context,  beamforming  can  be  used 
to  separate  signals  received  by  sensor  arrays  based  on  frequency,  wavenumber, 
and  velocity  (speed  as  well  as  direction)  of  propagation.  Both  the  transfer 
and  unit  impulse-response  functions  of  a  velocity  filter  are  two-dimensional 
functions  in  the  case  of  one-dimensional  arrays.  The  transfer  function 
involves  frequency  and  wavenumber  (due  to  spatial  sampling  by  equally  spaced 
sensors)  as  independent  variables,  whereas  the  unit  impulse  response  depends 
upon  time  and  location  within  the  array.  Two-dimensional  filtering  is  not 
limited  to  velocity  filtering  by  means  of  seismic  array.  Two-dimensional 
spatial  filters  are  frequently  used,  for  example,  in  the  interpretation  of 
gravity  and  magnetic  maps  to  differentiate  between  regional  and  local 
features.  Input  data  for  these  filters  may  be  observations  in  the  survey  of 
an  area  conducted  over  a  planar  grid  over  the  earth's  surface.  Two- 
dimensional  wavenumber  digital  filtering  principles  are  useful  for  this 
purpose.  Velocity  filtering  by  means  of  two-dimensional  arrays  may  be 
accomplished  by  properly  shaping  a  three-dimensional  response  function 
H(ki ,  k, .  (0)  .  Velocity  filtering  by  three-dimensional  arrays  may  be  accomplished 
through  a  four-dimensional  function  H(k^ ,  k2,  k3,  c*>)  as  explained  below. 
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Spatial  Arrays,  Beamf ormers,  and  FIR  Filters: 

A  propagating  plane  wave,  six,  t) ,  is,  in  general,  a  function  of  the 
three-dimensional  space  variables  (x,,x2,x3)  Ax  and  the  time  variable  t-  The 
4-D  Fourier  transform  of  the  stationary  signal  s(x,t)  is 


rm  f-  rm  ,  (1) 

S{&, <■>)=]  J  J  J  six.t)e  '  **  'dx3dx2dx3dt 

which  is  referred  to  as  the  wavenumber-frequency  spectrum  of  six,t),  and 
(k3 ,  k? .  k3)  A  k  denotes  the  wavenumber  variables  in  radians  per  unit  distance  and 
is  the  frequency  variable  in  radians  per  second.  If  c  denotes  the  velocity 
of  propagation  of  the  plane  wave,  the  following  constraint  must  be  satisfied 


+k!  = 


or 


If  the  4-D  Fourier  transform  of  the  unit  impulse  response  hix,  t)  of  a  4-D 
linear  shift-invariant  (LSI)  filter  is  denoted  by  Hik,  w)  ,  then  the  response 
y(x,t)  of  the  filter  to  s(x,t)  is  the  4-D  linear  convolution  of  hix,  t)  and 
six,  C) ,  which  is,  uniquely,  character ized  by  its  4-D  Fourier  transform 

Y(1 c,u)  =  Hik,  w)  Sik,  <*>)  •  ^ 


The  inverse  4-D  Fourier  transform,  which  forms  a  4-D  Fourier  transform  pair 
with  (1)  is 


1  rm  rm  rm  rm  kt*i\ 

six,  t)  *  i — I  iff  Sik,  w)  e  '  1  dk^dk^k^du) 

(2tt)  *  J— J— 


(3) 


It  is  noted  that  S(^,w)  in  (1)  is  product  separable,  i.e.  expressible 
in  the  form 


Slk.w)  *  S,  (/f1)S?  (/c2)Sj(A:j)S4(w)  , 


(4) 


where  each  function  on  the  right-hand  side  is  a  univariate  function  of  the 
respective  independent  variable,  if  and  only  if  six,  t)  in  (1)  is  also  product 
separable.  In  beamforming,  Siiki)  in  (4)  would  be  the  far-field  beampattern 
of  a  linear  array  along  the  Xj-axis.  For  example,  the  normalized  beampattern 
of  a  uniformly  weighted  (shaded)  linear  array  of  length  L  is 


Sik, 0) 


sin. 


JcLsinS 

2 


sin0 


(5) 


where  A,  =  |  — ^-j  is  the  wavelength  of  the  propagating  plane  wave  and  0  is  the 

angle  of  arrival  at  array  site  as  shown  in  Figure  1.  Note  that  0  is 
explicitly  admitted  as  a  variable  in  Sik.Q)  to  allow  for  the  possibility  that 
for  a  fixed  wavenumber,  the  beampattern  could  be  plotted  as  a  function  of  the 
angle  of  arrival.  In  that  case,  when  0  is  zero,  the  wave  impinges  the  array 
broadside  and  the  normalized  beampattern  evaluates  to  unity. 

The  counterpart,  in  aperture  and  sensor  array  processing,  of  the  use  of 


b> 
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window  functions  in  spectral  analysis  for  reduction  of  sidelobes  is  the  use  of 
aperture  shading.  In  aperture  shading,  one  simply  multiplies  a  uniformly 
weighted  aperture  by  the  shading  function.  The  resulting  beampattern  is, 
then,  simply  the  convolution  of  the  beampattern  of  the  uniformly  shaded 
volumetric  array  and  the  beampattern  of  the  shading  function.  Fourier 
transform  relationship  between  the  stationary  signal  s(ar,  t)  and  the 
wavenumber  frequency  spectrum  S(k,  o)  allows  one  to  exploit  high  resolution 
spectral  analysis  techniques  for  the  high  resolution  estimation  of  the 
direction  of  arrival  [Pillai,  1989] . 

Discrete  Arrays  for  Beamforming: 

An  array  of  sensors  could  be  distributed  at  distinct  points  in  space  in 
various  ways.  Line  arrays,  planar  arrays,  and  volumetric  arrays  could  be 
either  uniformly  spaced  or  nonuniformly  spaced,  including  the  possibility  of 
placing  sensors  randomly  according  to  some  probability  distribution  function. 
Uniform  spacing  along  each  of  the  coordinate  axis  permits  one  to  exploit  the 
well-developed  multidimensional  signal  processing  techniques  concerned  with 
filter  design,  DFT  computation  via  FFT  and  high  resolution  spectral  analysis 
of  sampled  signals  [Dudgeon,  1977],  Nonuniform  spacing,  sometimes,  might  be 
useful  for  reducing  the  number  of  sensors,  which,  otherwise,  might  be 
constrained  to  satisfy  a  maximum  spacing  between  uniformly  placed  sensors  to 
avoid  grating  lobes  due  to  aliasing,  as  explained  later.  A  discrete  array, 
uniformly  spaced,  is  convenient  for  the  synthesis  of  a  digital  filter  or 
beamformer  by  the  performing  of  digital  signal  processing  operations  (namely 
delay,  sum,  and  multiplication  or  weighting)  on  the  signal  received  by  a 
collection  of  sensors  distributed  in  space.  The  sequence  of  the  nature  of 
operations  dictate  the  types  of  beamformer.  Common  beamforming  systems  are  of 
the  straight  summation,  delay-and-sum,  and  weighted  delay-and-sum  types.  The 
geometrical  distribution  of  sensors  and  the  weights  wt  associated  with  each 
sensor  are  crucial  factors  in  the  shaping  of  the  filter  characteristics.  In 
the  case  of  a  linear  array  of  N  equispaced  sensors,  which  are  spaced  D  units 
apart,  starting  at  the  origin  x.,=0,  the  function 


w  =  jjZ  w*e 
lv  0 


Jk-_nD 


(6) 


becomes  the  array  pattern,  which  may  be  viewed  as  the  frequency  response 
function  for  a  finite  impulse  response  (FIR)  filter,  characterized  by  the  unit 
impulse  response  sequence  .  In  the  case  when  wn- 1,  (6)  simplifies  to 


(7) 


If  the  //  sensors  are  symmetrically  placed  on  both  sides  of  the  origin, 
including  one  at  the  origin,  and  the  sensor  weights  are  (v„  =  l,  then  the  linear 
array  pattern  becomes 


sin 


sin 


k.ND 

2 

kxD 

~z~ 


For  planar  arrays,  direct  generalizations  of  the  preceding  linear  array 
results  can  be  obtained.  To  wit,  if  the  sensors  with  unity  weights  are 

located  at  coordinates  (kD.iD)  ,  where  k  -  0 ,  ±  1 ,  ±  2 ,  .  .  .  ,  j ,  anci 
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{  =  0,  ±1,  ±2,  .  .  .  ,  for  odd  integer  values  of  AT  and  M  then  the  array 

pattern  function  becomes 

.  (t1)  m 

Kki‘ki)  =  E  E  exp{-j(k1lcD*lcJtZ}} 

*-W) 
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sinj 
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,  2  j 
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k2MD) 
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1  kl° 
l  2 
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K  2 

Routine  generalizations  to  3-D  spatial  arrays  are  also  possible.  The  array 
Pattern  functions  for  other  geometrical  distributions  may  also  be  routinely 
generated.  For  example,  if  unit  weight  sensors  are  located  at  the  six 
vertices  and  the  center  of  a  regular  hexagon  each  of  whose  sides  is  D  units 
long,  then  the  array  pattern  function  can  be  shown  to  be 

f^k1,k2)  -  y  1  +  2cos  k2D  *  4cos  cos  ~k~.  . 

The  array  pattern  function  reveals  how  selective  a  particular 
beamforming  system  is.  In  the  case  of  a  typical  array  function  shown  in  (7), 
the  beamwidth,  which  is  the  width  of  the  main  lobe  of  the  array  pattern,  is 
inversely  proportional  to  the  array  aperture.  Because  of  the  periodicity  of 

the  array  pattern  function,  the  main  lobe  is  repeated  at  intervals  of  . 

These  repetitive  lobes  are  called  grating  lobes,  whose  existence  may  be 
interpreted  in  terms  of  spatial  frequency  aliasing  resulting  from  a  sampling 
interval  D  due  to  the  N  receiving  sensors  located  at  discrete  points  in 
space.  If  the  spacing  D  between  sensors  satisfies 

(10) 


where  A  is  the  smallest  wavelength  component  in  the  signal  received  by  the 
array  of  sensors,  then  the  grating  lobes  have  no  effect  on  the  received 
signal.  A  plane  wave  of  unit  amplitude  which  is  incident  upon  the  array  at 
bearing  angle  0  degrees  as  shown  in  Figure  1  produces  outputs  at  the  sensors 
given  by  the  vector 

Sq  -  [exp ( j'O)  expfj^Dsine).  .  .  exp^/Cj  (W-l)Usin0)] c  (11) 


where  kx  -  is  the  wavenumber.  In  array  processing,  the  array  output  y# 

may  be  viewed  as  the  inner  product  of  an  array  weight  vector  jy  and  the 
steering  vector  .  Thus,  the  beamformer  response  along  a  direction 
characterized  by  the  angle  0  is 

N-l 

y#  =  (#'£*)  =  E  wkex p(jk1kDsir&) .  (12) 

#”•0 

The  Leamforming  system  is  said  to  be  robust  if  it  performs  satisfactorily 
despite  certain  perturbations  [Ahmed  and  Evans,  1982].  It  is  possible  for 
each  component  of  £,  to  belong  to  an  interval  » s»«  +$*ej  *  and  a  robust 

beamformer  will  require  the  existence  of  at  least  one  weight  vector  iy  which 


-4— 


will  guarantee  the  output  y#  to  belong  to  an  output  envelope  for  each  fig  m 
the  input  envelope.  The  robust  beamforming  problem  can  be  translated  into  an 
optimization  problem,  which  may  be  tackled  by  minimizing  the  value  of  the 
array  output  power 


P(6)  =  ac(0)Pisr(e) 


(13) 


when  the  response  to  a  unit  amplitude  plane  wave  incident  at  the  steering 
direction  0  is  constrained  to  be  unity,  i.e.  kc(0)s(0)  =1»  and  R  is  the 
additive  noise  corrupted  signal  autocorrelation  matrix.  The  solution  is 
called  the  minimum  variance  beamformer  and  is  ‘given  by 


(6)  -  S'1^) 

Xftv  '  sc(0)£-\s(0) 


(14) 


and  the  corresponding  power  output  is 

Pw(0)  = 


Sc(0)P-’fL(0) 


(15) 


The  minimum  variance  power  as  a  function  of  0  can  be  used  as  a  form  of  the 
data-adaptive  estimate  of  the  directional  power  spectrum.  However,  in  this 
mode  _f  solution,  the  coefficient  vector  is  unconstrained  except  at  the 
steering  direction.  Consequently,  a  signal  tends  to  be  regarded  as  an 
unwanted  interference  and  is,  therefore,  suppressed  in  the  beamformed  output 
unless  it  is  almost  exactly  aligned  with  the  steering  direction.  Therefore, 
it  is  desirable  to  broaden  the  signal  acceptance  angle  while  at  the  same  time 
preserving  the  optimum  beamformer' s  ability  to  reject  noise  and  interference 
outside  this  region  of  angles.  One  way  of  achieving  this  is  by  the 
application  of  the  principle  of  superdirectivity. 

Discrete  Arrays  and  Polynomials : 

It  is  common  practice  to  relate  discrete  arrays  to  polynomials  for  array 
synthesis  purposes  (Steinberg,  1976]  .  For  volumetric  equispaced  arrays  (it  is 
only  necessary  that  the  spacing  be  uniform  along  each  of  the  coordinate  axis 
so  that  the  spatial  sampling  periods  Dt  and  Dj  along,  respectively,  the  j ch 
and  jch  coordinate  axes  could  be  different  for  i  *  j)  the  weight  associated 
with  sensors  located  at  coordinate  (1,0, ,  i?D,,  i2D3)  is  denoted  by  w(i,,i2,i3)  . 
The  function  in  the  complex  variables  z, ,  z2 ,  ,  and  z3  that  is  associated  with 
the  sequence  {v(i, ,  i2,  i3)  }  is  the  generating  function  ft  the  sequence  and  is 
denoted  by 


fflz1,z2.zJ)  i2,  i-})  zi'zt'z^ .  (15) 

■h  •»> 

In  the  electrical  engineering  and  in  the  geophysics  literature,  the  generating 
function  W(z,,  z2,  z3)  is,  sometimes,  called  the  z-transform  of  the  sequence 
{w  (i, ,  i7,  i3) }  .  When  there  are  a  finite  number  of  sensors,  a  realistic 
assumption  for  any  physical  discrete  array,  Wlz,,z2,z3)  becomes  a  trivariate 
polynomial.  In  the  special  case  when  w(i1,iz,i3)  is  product  separable,  the 
polynomial  W(z3,z2,z3)  is  also  product  separable.  Particularly,  this 
separability  property  holds  when  the  shading  is  uniform,  i.e.  w (1, ,  i2,  i3)  =  1. 
When  the  support  of  the  uniform  shading  function  is  defined  by 

i,  =0,1 . A/,-1,  i2  -  0 , 1 . A/2-l,  and  i3  *  0,1 . A/3  -  1 ,  the  associated 

polynomial  becomes 


(17) 


If.'l  JT,-1  (1,-1  J 

W{z^.z2.zz)  =  £  £  X!  z*'z*  'z*  =  n 

il«0  a2»0  ij»0  i»l 


In  this  case,  all  results  developed  for  the  synthesis  of  linear  arrays  become 
directly  applicable  to  the  synthesis  of  volumetric  arrays.  For  a  linear 
uniform  discrete  array  composed  of  2V  sensors  with  intersensor  spacing  23, 
starting  at  the  origin  and  receiving  a  signal  at  a  known  fixed  wavenumber  k, 
at  a  receiving  angle  0,  the  far-field  beampattern 


S  ( k, ,  0 )  A  5(6) 


N-  1 


^k^cD ,sin6 


may  be  associated  with  a  polynomial  Zxr,  by  setting  z,  =  e-2*>c-Bln  .  This 

r«0 

polynomial  has  all  its  zeros  on  the  unit  circle  in  the  z, -plane.  If  the  array 
just  considered  is  not  uniform  but  has  a  weighting  factor  wr,  for 
r  =  0,1, ...,27,-1,  the  space  factor, 

0(0)  A  V  wreJk'D'caXn* 

r?0 


W,-l 

may  again  be  associated  with  a  polynomial  wrz,r  •  BY  the  pattern 

r-o 

multiplication  theorem,  it  is  possible  to  get  the  polynomial  associated  with 
the  total  beam  pattern  of  an  array  with  weighted  sensors  by  multiplying  the 
polynomials  associated  with  the  array  element  pattern  and  the  polynomial 
associated  with  the  space  factor  0(0)  .  The  array  factor  jo(0)  |2  may  also  be 
associated  with  the  polynomial  spectral  factor 


|O<0)  |2 


Wi-l 

-J 

r-0 


5T  WrZ^  ^2  Wr(ZlY > 


r»  0 


(18) 


where  the  weighting  (shading)  factor  is  allowed  to  be  complex.  Uniformly 
distributed  apertures  and  uniformly  spaced  volumetric  arrays  which  admit 
product  separable  sensor  weightings  can  be  treated  by  using  the  well-developed 
theory  of  linear  discrete  arrays  and  their  associated  polynomial.  When  the 
product  separability  property  does  not  hold,  scopes  exist  for  applying  results 
from  multidimensional  systems  theory  [Bose,  1982)  concerning  multivariate 
polynomials  to  the  synthesis  problem  of  volumetric  arrays. 

Velocity  Filtering: 

Combination  of  individual  sensor  outputs  in  a  more  sophisticated  way 
than  delay-and-sum  technique  leads  to  the  design  of  multichannel  velocity 
filters  for  linear,  planar,  as  well  as  spatial  arrays.  Consider,  first,  a 
linear  (1-D)  array  of  sensors,  which  will  be  used  to  implement  velocity 
discrimination.  The  pass  and  rejection  zones  are  defined  by  straight  lines  in 
the  (k, ,  caj-plane,  where 


* 


p 

V 


Ci) 

( v/sin0) 


is  the  wavenumber,  o  the  angular  frequency  in  radians/sec,  V  the  apparent 
velocity  on  the  earth's  surface  along  the  array  line,  v  is  the  velocity  of 
wave  propagation,  and  0  is  the  horizontal  arrival  direction.  The  transfer 
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function 


H((i) ,  k ,) 


otherwise, 


of  a  "pie-slice"  or  "fan"  velocity  filter  [Bose,  1985]  rejects  totally 

wavenumbers  outside  the  range,  -  s  kt  s.  -L^i,  and  passes  completely 

wavenumbers  defined  within  that  range.  Thus,  the  transfer  function  defines  a 
high-pass  filter  which  passes  signals  with  apparent  velocities  of  magnitude 
greater  than  v  at  a  fixed  frequency  w •  If  the  equi-spaced  sensors  are  D 
units  apart,  the  spatial  sampling  results  in  a  periodic  wavenumber  response 

with  period  k±  -  -uL  .  Therefore,  for  a  specified  apparent  velocity  V,  the 

resolvable  wavenumber  and  frequency  bands  are,  respectively,  - 

and  -  — —  s  w  £,  where  represents  the  folding  frequency  in  radians/sec. 

2D  2D  2D 


Linear  arrays  are  subject  to  the  limitation  that  the  source  is  required 
to  be  located  on  the  extended  line  of  sensors  so  that  plane  wavefronts 
approaching  the  array  site  at  a  particular  velocity  excite  the  individual 
sensors,  assumed  equi-spaced,  at  arrival  times  which  are  also  equi-spaced.  In 
seismology,  the  equi-spaced  interval  between  successive  sensor  arrival  times 

is  called  a  move-out  or  step-out  and  equals  ^s^-n—  =  — -  However,  when  the 

v  V 

sensor-to-source  azimuth  varies,  two  or  more  independent  signal  move-outs  may 
be  present.  Planar  (2-D)  arrays  are  then  required  to  discriminate  between 
velocities  as  well  as  azimuth.  Spatial  (3-D)  arrays  provide  additional  scope 
to  the  enhancement  of  discriminating  capabilities  when  sensor / source  locations 
are  arbitrary.  In  such  cases,  an  array  origin  is  chosen  and  the  mth  sensor 
location  is  denoted  by  a  vector  (xtm  x2m  x2m) c  and  the  frequency  wavenumber 
response  of  an  array  of  N  sensors  is  given  by 


i  w  3 

=  -±- £  Hm( o>)exp  T  -j2nkixil 

/v  m«  i  i-l 


where  denotes  the  frequency  response  of  a  filter  associated  with  the 

mth  recording  device  (sensor) .  The  sum  of  all  N  filters  provides  flat 
frequency  response  so  that  waveforms  arriving  from  the  estimated  directions  of 
arrival  at  estimated  velocities  are  passed  undistorted  and  other  waveforms  are 
suppressed.  In  the  planar  specialization,  the  2-D  array  of  sensors  lead  to 
the  theory  of  3-D  filtering  involving  a  transfer  function  in  the 
f requency/wavenumber  variables  f,ky  and  k?  .  The  basic  design  equations  for 
the  optimum,  in  the  least-mean-square  error  sense,  frequency  filters  have  been 
developed  [Burg,  1964].  This  procedure  of  Burg  can  be  routinely  generalized 
to  the  4-D  filtering  problem  mentioned  above. 
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Abstract 
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1  Introduction 


Sengupta  and  litis  [1]  formulated  the  computation  of  the  a  posteriori  probabilities  for  the 
joint  probabilistic  data  association  filter  (JPDAF)  [2]  as  a  constrained  minimization  prob¬ 
lem.  This  technique  based  on  the  use  of  neural  networks  was  also  extended  [3]  to  apply 
to  maneuvering  targets.  The  a  posteriori  probability  (3 j  (for  j  ^  0)  is  the  probability  that 
measurement  j  originated  from  target  t  and  (3lQ  is  the  probability  that  none  of  the  received 
measurements  originated  from  target  t.  By  comparison  with  the  traveling  salesman  problem 
(TSP),  they  proposed  a  Hopfield  neural  network  [4]  to  approximately  compute  the  /3'5s  and 
called  this  neural  network  probabilistic  data  association  (NPDA).  In  fact,  is  approximated 
by  the  output  voltage  V*  of  a  neuron  in  an  (m  +  1)  x  n  array  of  neurons,  where  m  is  the 
number  of  measurements  and  n  is  the  number  of  targets.  Sengupta  and  litis  claimed  that  the 
performance  of  the  NPDA  was  close  to  that  of  the  JPDAF  in  situations  where  the  numbers 
of  measurements  and  targets  were  in  the  ranges  of  3  to  20  and  2  to  6,  respectively.  The 
success  of  the  NPDA  in  their  examples  was  credited  to  the  accurate  emulation  of  all  the 
properties  of  the  JPDAF  by  the  NPDA. 

Since  there  are  no  strict  guidelines  for  choosing  the  constant  coefficients  of  the  energy 
function  for  the  Hopfield  neural  networks  [4],  these  coefficients  in  the  NPDA  were  selected 
essentially  by  trial-and-error  [1],  However,  with  the  chosen  coefficients  in  [1],  the  architecture 
of  the  neural  network  in  [1]  degenerates  into  individual  columns  with  identical  connections 
and  the  input  currents  throughout  the  neural  network  are  almost  of  the  same  amplitude. 
Furthermore,  the  neural  network  in  [1]  has  a  strong  tendency  to  converge  to  the  Vj's  which 
are  close  to  l/(m  +  1)  with  the  given  initial  values  for  the  V-'s.  In  [3],  a  least-squares 
approach  was  developed  to  optimize  the  last  two  coefficients  of  the  energy  function.  The 
authors  pointed  out  that  this  local  optimization  of  the  two  coefficients  is  quite  sensitive  to 
ra,  the  number  of  targets.  They  chose  the  two  coefficients  by  averaging  their  values  for  the 
n  =  2,  n  =  4,  and  n  =  6  cases.  As  a  result,  the  coefficients  of  the  energy  function  in  [3]  are 
slightly  different  from  the  ones  in  the  energy  function  of  [1].  Because  of  the  similarity  between 
the  energy  functions  in  [1,  3] ,  the  discussion  in  this  correspondence  will  be  focused  on  the 


2 


energy  function  in  [1].  Similar  discussions  on  the  energy  function  in  [3]  are  documented  in 
detail  in  [5]. 

In  Section  2,  the  Hopfield  neural  network  used  in  [1]  is  briefly  reviewed  and  some 
comments  are  made  on  the  assumptions  which  were  used  to  set  up  the  energy  function  in  [1]. 
In  Section  3,  some  criticisms  on  the  implementation  of  the  neural  network  in  [1]  are  given. 
Conclusions  are  summarized  in  Section  4. 


2  Review  of  the  Energy  Function  in  the  NPDA  and  Comments 

Suppose  there  are  n  targets  and  m  measurements.  The  energy  function  used  in  [1]  is  repro¬ 
duced  below. 

A  m  n  n  n  n  m  m  p  n  m 

w  =  ^EEE^r  +  TlEEw4E(E^-i)!  + 


2  z- 

*  j= 0 4=1 r  =  1 
T*t 


i=l>=0/  =  0 

I  £  i 


t= 1  j=0 


r\  m  n  jp  m  n  n  m 

t  £  EW  -  d)s  +  #EE£(»7-E  rf  )2- 


2  4 

z  j=Ol=l 


(1) 


j=Ot=l  r  =  1 
r  /  t 


l  =  0 


In  (1),  Vf  is  the  output  voltage  of  a  neuron  in  an  (m  +  1)  x  n  array  of  neurons  and  is 
the  approximation  to  the  a  posteriori  probability  /3‘  in  the  JPDAF  [2].  This  a  posteriori 
probability,  in  the  special  case  of  the  PDAF  [2]  when  the  probability  PG  that  the  correct 
measurement  falls  inside  the  validation  gate  is  unity,  is  denoted  by  p'.  Actually,  PG  is  very 
close  to  unity  when  the  validation  gate  size  is  adequate,  as  shown  in  Table  5-1  in  [2].  In  (1), 
A,  £?,  C ,  D,  and  E  are  constants. 

In  the  NPDA,  the  connection  strength  matrix  is  a  symmetric  matrix  of  order  n(m+l), 
With  the  given  energy  function  Eoap  in  (1),  the  connection  strength  T*{  from  the  neuron 
at  location  (r,  /)  to  the  neuron  at  location  (t,j)  is 

—(C  +  D  +  E(n  —  1)),  if  t  =  r  and  j  =  l  “self  feedback”, 

—  A,  lit  ^  t  and  j  =  l  “row  connection”, 

— (B  +  C),  lit  —  r  and  j  ^  /  “column  connection”, 

0,  lit  ^  r  and  j  ^  l  “global  connection”. 


Tjl  =  i 


(2) 
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The  input  current  /j  to  the  neuron  at  location  (f,j>),  for  t  =  1,2, . . .  ,ra,  and  j  =  0,1,..., m, 
is 

/'  =  C  +  (.D  +  E)p‘i  +  B(n -l -■£,;).  (3) 

r=l 

Clearly  from  (2)  and  (3),  the  input  current  Ij  but  not  the  connection  strength  T‘[  depends 
on  the  p1-' s,  which  are  computed  from  the  measurements  that  comprise  the  input  data. 
Ironically,  in  the  neural  network  for  the  TSP  [4],  only  the  connection  strengths  depend  on 
the  input  data  which,  in  this  case,  are  the  distances  between  pairs  of  cities. 

In  order  to  justify  the  first  two  terms  of  Edap  in  (1),  the  authors  of  [1]  claimed  that 
the  dual  assumptions  of  no  two  returns  from  the  same  target  and  no  single  return  from  two 
targets  are  consistent  with  the  presence  of  a  dominating  Vj  in  each  row  and  each  column 
of  the  (m  +  1)  x  n  array  of  neurons.  However,  these  assumptions  are  not  constraints  on 
the  values  of  the  /3j’ s  in  the  original  JPDAF.  Those  assumptions  should  be  used  only  in  the 
generation  of  the  feasible  data  association  hypotheses,  as  pointed  out  in  [2,  p.  224].  As  a 
matter  of  fact,  there  could  be  two  /3j’s  of  comparable  magnitude  in  the  same  row  and  in  the 
same  column  as  shown  in  chapter  4  of  [5].  Therefore,  the  presence  of  a  dominating  Vj  in 
each  row  and  each  column  is  not  a  property  of  the  JPDAF. 

The  third  term  of  Edap  is  used  to  constrain  the  sum  of  the  Vj' s  in  each  column  to 

m  m 

unity  i.e.  )  Vj  =  1.  This  constraint  is  consistent  with  the  requirement  that  =  1  in 

j=o  j= o 

both  the  JPDAF  and  the  PDAF  [2].  Therefore,  this  constraint,  by  itself,  does  not  permit  us 

to  infer  whether  the  /3j's  are  from  the  JPDAF,  or  from  the  PDAF.  The  assumption  used  to 

set  up  the  fourth  term  is  that  this  term  is  small  only  if  Vj  is  close  to  p‘J}  in  which  case  the 

neural  network  simulates  more  closely  the  PDAF  for  each  target  rather  than  the  intended 

JPDAF  in  the  multitarget  scenario.  Finally,  the  fifth  term  is  supposed  to  be  minimized  if  Vj 

is  not  large  unless  for  each  r  /  f  there  is  a  unique  l  j  such  that  p [  is  large.  Unfortunately, 

this  constrained  minimization  may  not  be  possible  as  shown  in  [5].  This  is  consistent  with 

the  heuristic  nature  of  the  derivation  of  the  energy  function  in  [1],  which  could  lead  to  the 

problems  in  the  implementation  of  the  NPDA  as  discussed  next. 
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3  Criticism  of  the  Implementation  of  the  NPDA 

Since  there  are  no  guidelines  for  choosing  A,  B,  C,  D,  and  E  in  (1),  these  coefficients  are 
usually  set  through  simulation.  In  the  implementation  of  the  NPDA  [1],  the  simulations 
provided  the  following  set  of  values  for  A,  B ,  C,  D,  and  E  for  tracking  two  to  six  targets 
with  three  to  twenty  measurements. 

A  =  0  jB  =  40  C  =  1000 
D  =  30  E  =  10 

After  substituting  the  above  set  of  values  into  (2),  the  only  connections  between  neurons 
which  are  nonzero  are  those  due  to  “self  feedback”  and  “column  connection”,  as  evident  from 
(4)  below. 

—(1030  +  10(n  —  1)),  if  t  =  r  and  j  =  l  “self  feedback” 

0,  if  t  ^  r  and  j  =  l  “row  connection” 

(4) 

—  1040,  if  t  =  r  and  j  ^  l  “column  connection” 

0,  if  t  ^  t  and  j  /  l  “global  connection”. 

Since  0  <  P'  <  1,  the  upper  and  lower  bounds  of  the  input  current  /j  in  (3)  can  be 
obtained. 

1000  <  I)  <  1030  +  10(n  -  1).  (5) 

The  upper  bounds  of  the  input  current  /j  in  [1]  are  1040,  1060,  and  1080  for  n  =  2,4,6, 
respectively.  Clearly,  the  maximum  range  of  variation  of  the  input  current  is  only  a 
small  percentage  of  the  value  of  C .  Therefore,  the  input  current  /j  is  dominated  by  C  for 
n  =  2,4,6.  In  other  words,  the  p^'s  do  not  have  much  influence  on  the  solution  obtained 
from  the  neural  network. 

One  may  argue  that  C  does  not  dominate  the  behavior  of  the  neural  networks  if  the 
dynamic  equation  (32)  in  [l]  is  considered.  For  the  purpose  of  the  discussion,  this  equation 
with  A  =  0  is  reproduced  below. 


5 


(6) 


dul 

j 

da 


ul.  m  /  m  \ 

~J--B'EV?-C[J2V?-l)-(D  +  E{n- l))V‘ 

s°  1  =  0  \i= o  / 


+(D  +  E)P'j  +  E(n-  !-£>;). 


T  =  1 


In  (6),  so  is  the  time  constant  of  the  neural  network  and  u‘  is  the  input  voltage  of  the  neuron 
located  at  (t,j)  in  the  array  of  neurons.  When  the  neural  network  is  initialized  so  that  each 
neuron  approximately  outputs 


v  =  ,J_ 

3  m  +  1 


(7) 


rn 

the  third  term  is  almost  zero  because  ^  Vf  is  close  to  unity.  However,  we  must  keep  in  mind 

1=0 

that  this  neural  network  has  not  been  developed  for  implementation  on  a  digital  computer. 
Instead,  it  has  been  developed  for  possible  implementation  using  an  analog  network.  If  it  is 
implemented  on  a  digital  computer,  (6)  is  replaced  by  (38)  in  [1],  which  is  reproduced  below 
for  ready  reference. 


“it* + 1)  =  f;  v<‘ -  fix  - 'l  -oD+£(»-i))v? 

Sq  1  * 


l  =  0 


w=o 


+<(*> +  £>}  + C£(«  - 1  -  £>;), 


(8) 


T=  1 


In  (8),  (  is  the  iteration  step  and  12  multiplications  are  sufficient  for  implementing  each 
iteration.  A  total  of  400  iterations  was  used  to  guarantee  the  convergence  of  the  neural 
network  to  a  satisfactory  solution  [1].  There  are  altogether  n(m  +1)  equations  like  (8)  for 
tracking  n  targets  with  m  measurements.  Therefore,  the  total  number  of  multiplications 
required  to  compute  the  /3j’s  does  not  exceed  n(m  +  1)  *  12  *  400.  Comparing  this  number 
with  the  entries  in  Tables  2.4,  3.3,  and  3.4  in  [5],  we  can  conclude  that  the  implementation 
of  the  neural  network  on  a  computer  is  far  more  computationally  expensive  than  the  direct 
computation  of  the  /3j’s.  If  this  neural  network  is  implemented  on  an  analog  circuit,  the 
dynamic  equation  of  the  network  is 


u\ 


30  r= 1 1=0 


(9) 
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as  shown  in  (3)  of  [4],  Therefore,  in  this  example  because  of  the  7j  term  in  (9),  it  is  seen 
from  (5)  that  C  becomes  a  dominan,  .  ctor  in  the  neural  network.  Since  the  range  of 
variation  of  Jj  is  only  a  small  percentage  of  C ,  any  perturbation  in  the  implementation  of 
the  current  source  due  to  noise  would  contribute  to  the  output  voltage  Vj ,  which  is  supposed 
to  approximate  /9j.  Since  C  is  significantly  larger  than  A ,  B,  D,  and  E,  the  network  has  a 

m 

strong  tendency  to  be  locked  into  a  state  where  V,1  =  1.  With  the  above  initialization, 

/=o. 

identical  connections,  and  nearly  uniform  input  currents,  it  is  even  easier  for  the  network  to 
converge  to  the  normalized  V-'s  which  are  close  to  1  /(m  +  1).  Therefore,  it  is  possible  for 
some  Vj's  to  be  not  0  even  if  the  /3j’s  are  0.  As  a  result,  the  state  of  target  t  is  updated  in 
part  by  measurement  j  which  is  not  inside  the  validation  gate  of  target  t. 

Finally,  it  is  worth  mentioning  that  the  example  given  in  [1]  for  tracking  six  targets  is 
not  actually  a  six-target-tracking  problem  from  the  JPDAF  point  of  view.  In  multitarget 
tracking,  targets  are  grouped  into  clusters  before  applying  either  the  JPDAF  or  the  PDAF. 
Targets  are  in  the  same  cluster  if  there  is  at  least  one  measurement  in  each  of  the  intersections 
of  the  validation  gates  of  these  targets.  The  JPDAF  is  only  applied  to  a  cluster  which  contains 
more  than  one  target.  Otherwise,  the  PDAF  is  applied.  If  the  targets  in  the  examples  given 
in  [1]  are  grouped  into  clusters,  the  largest  size  of  a  cluster  of  targets  is  3.  In  such  cases,  the 
s  can  be  easily  computed  by  the  algorithms  proposed  in  chapters  2  and  3  of  [5]  and  also 
reported  in  [6]. 

4  Summary 

It  might  be  a  good  idea  to  use  the  neural  network  to  approximately  compute  the  a  posteriori 
probability  /3j  in  the  JPDAF.  However,  the  neural  network  developed  in  [1]  has  been  shown 
to  have  improper  energy  functions.  This  resulted  from  misinterpretations  of  the  properties 
of  the  JPDAF  which  the  network  was  supposed  to  emulate.  Furthermore,  improper  choices 
of  the  constant  coefficients  in  the  energy  function  in  [1]  make  the  situation  worse.  Therefore, 
further  study  of  the  properties  of  the  JPDAF  is  needed  to  construct  a  better  energy  function 
for  the  NPDA. 
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Abstract 

In  this  paper,  three  fast  algorithms  have  been  developed  to  strive  the  problem 
of  data  association  in  multitarget  tracking  in  clutter.  In  the  first  algorithm,  the 
problem  of  data  association  is  identified  as  an  exhaustive  search  problem  in  general. 
Subsequently,  a  mathematical  model  is  proposed  for  the  problem  of  data  associa¬ 
tion  in  the  joint  probabilistic  data  association  filter  (JPDAF).  Based  on  the  model, 
a  depth-first  search  (DFS)  approach  is  developed  for  the  fast  generation  of  data 
association  hypotheses  and  the  computation  of  the  conditional  probabilities  of  the 
hypotheses  in  the  JPDAF.  With  proper  preprocessing,  the  number  of  multiplications 
in  the  computation  of  the  conditional  probabilities  is  less  than  twice  the  number  of 
data  association  hypotheses.  Therefore,  the  computational  complexity  of  the  algo¬ 
rithm  is  analyzed  in  terms  of  the  number  of  data  association  hypotheses  in  the  worst 
case  and  also  in  the  average  case.  When  the  density  of  targets  is  moderate,  an  even 
more  efficient  algorithm  is  developed  to  directly  compute  a  posteriori  probabilities 
in  the  JPDAF  without  generating  the  data  association  hypotheses.  In  the  third  al¬ 
gorithm,  the  interference  among  closely  spaced  targets  is  simplified.  This  lead  us  to 
develop  an  approach  to  approximately  compute  the  a  posteriori  probabilities  in  the 
JPDAF. 

'This  research  has  been  sponsored  by  SDIO/IST  and  managed  by  the  Office  of  Naval  Research  under 
contract  NOOOH-86-K-0542 
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1  Introduction 


In  multitarget  tracking  in  clutter,  often  there  arc  more  than  one  measurement  available 
for  updating  the  state  of  a  single  target  .  To  solve  the  problem  of  data  association  between 
targets  and  measurements,  two  typical  approaches  have  been  reported  in  the  literature 
in  the  70’s.  One  is  called  the  track-oriented  approach  in  which  each  measurement  is 
assumed  to  have  originated  from  either  a  known  target  or  clutter,  as  in  Probabilistic 
Data  Association  Filter  (PDAF)  j  1,2]  and  Joint  Probabilistic  Data  Association  Filter 
(JPDAF)  1 1,31.  The  other  is  called  a  measurement-oriented  approach  in  which  each 
measurement  is  hypothesized  to  have  originated  from  either  a  known  target,  a  new 
target,  or  clutter  |4i.  In  both  approaches,  the  number  of  data  association  hypotheses 
could  increase  rapidly  with  the  increase  of  the  number  of  targets  and  the  number  of 
measurements.  Therefore,  the  computational  cost  in  generating  the  data  association 
hypotheses  would  overwhelmingly  dominate  in  a  multitarget  tracking  algorithm  when 
the  number  of  targets  and  the  number  of  measurements  are  relatively  large.  In  the 
PDAF,  the  computational  cost  for  data  association  is  reduced  drastically  by  isolating 
targets  from  each  other.  Only  the  measurements  which  lie  in  the  validation  gate  of  a 
target  are  considered  in  data  association.  A  validation  gate  is  a  specific  region  around 
the  predicted  position  of  a  target.  Since  the  PDAF  ignores  the  interference  from  other 
targets,  it  may,  sometimes,  cause  inistracking  of  closely  spaced  targets,  as  discussed 
by  Fortmann,  Bar-Shalom  and  Scheffe  1 3] .  This  difficulty  is  greatly  reduced  by  using 
the  JPDAF  [3j,  which  accounts  for  the  fact  that  a  measurement  which  falls  inside  the 
intersection  of  the  validation  gates  of  several  targets  could  have  originated  from  any 
one  of  these  targets  or  from  clutter.  In  the  JPDAF,  targets  are  divided  into  clusters. 
The  targets  are  in  the  same  cluster  if  there  is  at  least  one  measurement  inside  each  of 
the  intersections  of  their  validation  gates.  The  computational  cost  for  data  association 
is  reduced  by  grouping  targets  into  clusters.  The  number  of  different  data  association 
hypotheses  in  each  cluster  is  still  an  exponential  function  of  the  number  of  targets  in  the 
cluster.  Due  to  the  complexity  of  the  problem  of  data  association,  the  implementation 
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of  a  multitarget  tracking  algorithm  has  only  been  carried  out.  in  a  2-target  case  3,4,. 
except  in  5;.  In  j 5 i ,  an  initial  version  of  the  depth-first  search  (DFS)  algorithm  proposed 
in  this  paper  was  used  to  track  11  non-maneuvering  and  maneuvering  targets  through 
the  implementation  of  JPDA  filtering  and  smoothing.  Details  of  the  proposed  algorithm 
were  not  discussed  previously. 

To  reduce  significantly  the  computational  cost  for  data  association,  several  approx¬ 
imations  of  the  JPDAF  and  the  algorithm  in  |4j  have  been  reported  in  literature.  A 
first  version  of  the  approximate,  or  '’cheap'’,  JPDAF  was  published  in  j6j,  where  the 
probability  of  association  of  target  t  with  measurement  j  was  computed  by  an  ad  hoc 
formula.  The  author  of  G  developed  his  approach  further  in  i 71.  In  i7i,  the  performance 
of  various  versions  of  the  approximate  JPDAF’s  were  tested.  The  "cheap  JPDAF”  in 
G  performed  fairly  well  in  case  of  two  targets  but  failed  to  track  fo\ir  targets  [8].  In 
8  .  the  JPDAF  was  emulated  by  a  neural  network  which  is  capable  of  handling  two  to 
six  targets  and  three  to  twenty  measurements  with  a  set  of  carefully  chosen  coefficients 
for  the  network.  In  the  neural  network  approach,  the  a  posteriori  probabilities,  /Jj’s,  of 
association  of  target  t  with  measurement  j  were  computed  in  parallel  subject  to  the  as¬ 
sumption  that  one  out  of  all  3*  ’s  for  any  particular  target  t  is  significantly  larger  than  the 
remaining  j3*' s  for  the  same  target.  With  this  assumption,  the  resulting  approximation  of 
JPDAF  is  transformed  more  or  less  into  the  nearest-neighbor  JPDAF  7'.  Furthermore, 
in  the  neural  network  approach,  targets  are  not  grouped  into  clusters  before  calculating 
31  as  proposed  in  the  original  JPDAF  3j.  Alternatively,  Xagarajan.  e.t  al.  ,9j  arranged 
the  hypotheses  in  the  measurement-oriented  approach  |4i  in  a  special  order  so  that  the 
probabilities  of  the  hypotheses  are  proportional  to  the  product  of  certain  probability  fac¬ 
tors  already  evaluated.  The  algorithm  locates  the  N  globally  best  hypotheses  without 
evaluating  all  of  them.  For  tracking  multiple  maneuvering  targets  in  clutter,  the  authors 
of  |8j  introduced  the  approach  of  joint  probabilistic  data  and  maneuver  association  in 
lOi. 


In  this  paper,  three  fast  algorithms  have  been  developed.  In  the  first  algorithm, 
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a  mathematical  model  is  proposed  for  the  problem  of  data  association  in  general.  This 
model  is  based  on  exhaustive  search.  Specifically,  this  model  is  tailored  to  the  problem  of 
data  association  in  the  JPDAF.  Each  feasible  data  association  hypothesis  in  the  JPDAF 
is  equivalent  to  a  solution  to  the  exhaustive  search  problem  for  data  association.  Based 
on  the  proposed  model,  a  depth-first  search  (DFS)  algorithm  is  developed  to  generate 
data  association  hypotheses  efficiently  for  the  JPDAF.  The  manner  in  which  the  data 
association  hypotheses  are  generated  by  the  DFS  algorithm  suggests  an  approach  to 
further  reduce  the  number  of  multiplications  in  the  computation  of  the  conditional  prob¬ 
abilities  of  the  data  association  hypotheses.  In  fact,  after  proper  preprocessing,  the  total 
number  of  multiplications  in  the  computation  of  the  conditional  probabilities  for  the  data 
association  hypotheses  in  the  JPDAF  has  been  found  to  be  less  than  twice  the  number 
of  data  association  hypotheses  generated  by  exhaustive  search.  The  computational  com¬ 
plexity  of  the  proposed  DFS  algorithm  is  evaluated  in  terms  of  the  number  of  hypotheses 
both  in  the  worst  and  the  average  cases.  The  proposed  DFS  algorithm  is  most  suitable 
for  ground-based  surveillance  system  with  a  large  centralized  computational  capability, 
where  accurate  result  is  desired. 

However,  in  a  small  tracking  system  or  a  missile  guidance  system  where  computa¬ 
tional  capability  is  limited,  an  even  more  faster  algorithm  is  desired.  The  second  and 
the  third  algorithms  proposed  in  this  paper  are  designated  to  be  implemented  in  a  small 
system  or  a  system  with  multiprocessor.  In  the  second  algorithm,  the  a  posteriori  prob¬ 
abilities.  Jj’s.  are  computed  directly  without  generating  the  data  association  hypotheses 
in  the  JPDAF  when  either  the  density  of  targets  is  not  very  high,  or  the  size  of  the 
largest  cluster  of  targets  is  less  than  or  equal  to  4.  For  the  case  in  which  the  size  of 
cluster  of  targets  could  be  larger  than  4,  an  approximation  of  the  second  algorithm  is 
developed.  In  the  approximation,  only  the  interference  from  the  neighboring  targets  of 
target  t  are  considered  in  the  computation  of  /?j.  A  target  is  a  neighbor  of  target  t  if 
there  is  at  least  a  measurement  inside  the  intersection  of  the  validation  gates  of  the  two 
targets.  The  second  and  the  third  algorithms  are  equivalent  when  the  size  of  a  cluster 
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of  tarots  is  2.  Tlit'  important  feature  of  the  third  algorithm  is  that  3l  for  <'acli  target  is 
consulted  independently  as  in  the  TDAF  algoritlim. 


A  description  of  the  problem  of  interest  is  provided  in  Section  2.  In  this  section,  the 
.1PDAF  is  briefly  reviewed  to  explain  the  motivation  behind  this  research.  Furthermore, 
an  alternate  formula  for  the  conditional  probabilities  of  the  data  association  hypotheses 
is  derived  in  order  to  simplify  the  notation  in  the  development  of  the  three  algorithms. 
In  Section  3.1.  a  comparison  between  a  general  exhaustive  search  problem  and  the  one 
of  data  association  in  multitarget  tracking  is  discussed.  This  leads  to  a  mathematical 
model  for  the  problem  of  data  association.  In  Section  3.2  some  unique  features  of  the 
problem  of  data  association  are  pointed  out  as  a  prelude  to  a  specialized  DFS  algorithm 
for  generating  all  data  association  hypotheses.  In  Section  3.3.  the  advantage  of  the  DFS 
algorithm  is  further  demonstrated  in  the  computation  of  the  conditional  probabilities 
of  the  data  association  hypotheses.  I.i  Section  3.4.  the  computational  complexity  of 
the  proposed  DFS  algorithm  is  analyzed.  In  Section  4,  an  algorithm  is  proposed  to 
compute  directly  when  the  density  of  targets  is  not  very  high.  An  approximation 
of  the  direct  computation  algorithm  is  presented  in  Section  5.  Finally,  in  Section  6. 
computer  simulation  is  carried  out  for  tracking  multitarget  in  different  scenario. 


I 

I 

I 

\ 

I 
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2  Problem  Formulation 


Lot  in  and  u  bo  the  numbers  of  measurements  and  targets,  respectively,  in  a  particular 
cluster.  Jq  is  the  a  posteriori  probability  for  no  measurement  to  originate  from  target  t 
and  J‘  (j  r-  0)  is  the  a  posteriori  probability  for  measurement  j  to  originate  from  target 
t.  In  the  JPDAF.  the  computation  of  dj’s  begins  with  the  construction  of  a  so-called 
validation  matrix  11  for  the  n  targets  and  rn  measurements.  The  validation  matrix  fl  is 
a  m  x  (n  •+-  l)  rectangle  matrix  defined  as  [3] 

/ 

0  1  2  •  •  •  n 

•  .  .  .  \  ,  I 

1  -*-11  1 2  1  n  t 

1  -^‘21  --'22  ‘  ‘  -  X.’2n  2 

^  .  .  .  .  (  J'  l1) 

^  1  ^ml  ^m2  *  *  *  ^ mn  j 

where  1,  U/’j{  — -  I  if  measurement  j  is  inside  the  validation  gate  of  target  t  and 

jjp  =  0  if  measurement  j  is  outside  the  validation  gate  of  target  t  for  j  =  l,2,...,m, 

and  t  —  1.2 . n.  Based  on  the  validation  matrix  fl,  data  association  hypotheses  (or 

feasible  events  [31)  are  generated  subject  to  the  following  two  restrictions: 

1.  each  measurement  can  have  only  one  origin  (either  a  specific  target  or  clutter). 

2.  no  more  than  one  measurement  originates  from  a  target. 

This  leads  to  a  combinatorial  problem  where  the  number  of  data  association  hypotheses 
increases  exponentially  with  the  number  of  targets  and  the  number  of  measurements. 

Each  feasible  event  £  is  represented  by  a  hypothesis  matrix  J1  in  |3l.  fl  has  the  same 
size  as  the  validation  matrix  fl.  A  typical  element  in  fl  is  denoted  by  dj(  where  =  1 
only  if  measurement  j  is  hypothesized  to  be  associated  with  clutter  (/  0)  or  target 

t  (t  i  0).  After  each  fl  is  obtained,  the  conditional  probability  of  the  corresponding 
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data  association  hypothesis  or  feasible  event  is  calculated  by  a  formula  given  in  ,3  .  A 
simplified  version  of  this  formula  is  given  in  (2).  where  time  index  k  is  omitted. 


P{£(h)\Z)  -  Jj  pt  (2) 

6'  ah,  i 

for  j  1,2 . maud/  1.2....,n. 

where  Z  is  the  set  of  all  measurements  received  up  to  current  time  index  fc,  c  is  a 
normalizing  constant,  rna  is  the  number  of  targets  detected  in  this  feasible  event  <f, 
ujji  =  1  indicates  that  measurement  j  is  associated  with  target  t  in  the  event,  and 

pl  _  (  Nl7.‘-,0,S‘)PD  \lv,t  -=  1 

I  0  otherwise 

---  A(1  -PD)  =  P0.  (4) 

for  j  1, 2, . . . .  m  and  t  -  1.2 . n. 

In  (3)  and  (4),  A  is  the  clutter  density,  Pn  is  the  probability  of  detection,  and  N( z‘;0,  5') 
is  a  normal  distribution  density  function  with  zero  mean  and  covariance  matrix  S‘.  It 
is  understood  that  the  normalizing  constant  c  in  (2)  is  obtained  by  the  summation, 
P(£(n)|Z).  Therefore,  c  is  omitted  hereafter.  The  a  posteriori  probability/^'  is 
computed  from  the  conditional  probabilities  in  (2)  by 

,l‘  -  ■£  P[£(h)\Z)i,„  (5) 

fin, 

m 

4  --  i-Y.4  W 

for  j  —  1,2,  and!  =  1,2 _ ,n. 

In  order  to  have  a  better  understanding  of  how  each  /j'  is  calculated,  consider  an 
example.  Suppose,  there  are  two  targets.  In  one  radar  scan,  three  measurements  are 
received  and  suppose  one  of  these  falls  inside  the  intersection  of  the  validation  gates  of 
the  two  targets.  The  mathematical  representation  of  this  situation  may  be  described  by 
a  validation  matrix  12: 
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With  the  two  restrictions  given 
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(7) 


The  conditional  probability  P[£  (ttt)\Z)  can  be  easily  calculated  by  using  (2).  Subse¬ 
quently.  all  i3y s  may  be  computed.  For  example. 


W(n,)|z)  +  p(£(n2)|Z)H- W(n3)iz) 


The  problems  of  interest  here  in  the  implementation  of  the  JPDAF  for  any  number 
of  targets  and  measurements  are  the  efficient  generation  of  hypothesis  matrices  fi (  and 
fast  computation  of  j3‘ . 
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3  The  Depth-First  Search  (DFS)  algorithm 


3.1  Mathematical  Model  for  Data  Association 

As  pointed  out  earlier,  data  association  in  multitarget  tracking  is  a  combinatorial  prob¬ 
lem.  The  computational  cost  of  data  association  increases  exponentially  with  n,  the 
number  of  targets,  and  rn.  the  number  of  measurements.  The  efficiency  of  the  algorithm 
used  in  the  generation  of  the  data  association  hypotheses  is  especially  important  when 
n  and  rn  are  relatively  large.  In  order  to  develop  an  efficient  algorithm  to  generate  all 
data  association  hypotheses,  it  is  necessary  that  its  nature  be  well  understood.  So  a 
mathematical  model  is  developed  for  data  association. 

A  well-known  model  for  a  combinatorial  problem  is  called  exhaustive,  search  with 
certain  constraints [  1 1  j .  The  general  description  of  a  typical  exhaustive  search  problem 
is  following: 

There  are  rn  variables  Xj  ( j  -  1,2, ... ,  m).  The  value  of  each  X}  belongs 
to  a  set  Zj.  where  Z}  is  finite  and  linearly  ordered.  A  candidate  solution  for 
this  problem  is  a  m-tuple 

(X] .  A' 2, - A'm). 

The  objective  here  is  to  find  either  one  solution  or  all  solutions  each  of  which 
satisfies  the  imposed  constraints. 

An  efficient  algorithm  to  solve  the  above  problem  is  a  depth- first  search  (DFS)  procedure 
which  involves  backtracking  j  1 1 ) .  In  the  DFS  procedure,  a  solution  is  found  by  checking 
the  m  variables  in  a  m-tuple  one  by  one  from  left  to  right,  as  shown  in  Fig.  1.  In  Fig.  1, 
GET  X EXT(Z),  X})  is  a  logical  function.  It  is  true  only  if  X}  €  Z;,  X}  has  not  been 
used  in  the  DFS  procedure  yet.  and  X(,  X?,  . ..,  and  X}  satisfy  the  given  constraints. 
Otherwise,  it  is  false. 
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In  the  context  of  tracking  multitarget  in  clutter,  however,  data  association  can 
he  modelled  as  an  exhaustive  search  problem  with  a  set  of  proper  notations.  Let  X '} 

(j  1.2 . m)  denote  measurement  j.  The  value  of  Xj  identifies  the  target  or  clutter 

which  is  hypothesized  to  he  associated  with  measurement  j .  For  example,  X}  -  3  implies 
that  measurement  j  is  hypothesized  to  be  associated  with  target  3.  However,  if  A';  —  0, 
measurement  j  is  hypothesized  to  he  associated  with  clutter.  Z:  is,  therefore,  a  set 
of  clutter  denoted  by  0  and  those  targets  for  which  x]t  —  1  in  the  validation  matrix. 
Equivalently, 

Z-,  -  {t\<x]t  -  1},  j  —  1,2  and  t  -  0, 1,2,  ...,n.  (8) 

Note  that  since  u>j0  is  equal  to  1,  therefore  0  (E  Zj  for  j  -  1,2.. . .  ,  m.  If  measurement 
j  falls  inside  the  validation  gate  of  target  t,  t  6  Z}.  For  the  example  discussed  in  the 
previous  section,  the  Z} ’s  (j  -  1.2.3)  are: 

Z\  -  {0, 1};  Z2  =  {0.1.2};  Z,  -  {0,2}. 

In  the  JPDAF  scenario,  the  two  constraints  which  have  to  be  satisfied  for  a  feasible 
event  can  be  easily  translated  into  the  language  of  exhaustive  search  problem  for  data 
association.  Hence,  a  m-tuple 

( A  1 ,  A*2 ,  .  .  .  .  Xp,  ....  Xq ,  .  .  .  .  Xm  ) 
is  a  solution  if  the  following  two  constraints  are  satisfied. 

1.  If  p  f-  <j,  Xp  y  0,  and  Xq  /  0,  then  Xp  /-  A"7. 

2.  If  A'p  -  Xq  and  p  ^  q,  then  Xv  —  Xq  —  0. 

All  data  association  hypotheses  can  be  generated  by  solving  the  exhaustive  search 
problem  described  above.  In  the  next  section,  a  DFS  algorithm  will  be  discussed  in  detail 
for  generating  all  data  association  hypotheses. 
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Figure  1  The  flowchart  of  DFS  procedure 
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3.2  Generation  of  Data  Association  Hypotheses 

Ill  the  previous  subsection,  the  similarities  between  the  problem  of  exhaustive  search  and 
that  of  data  association  have  been  discussed  and,  subsequently,  a  mathematical  model 
lias  been  proposed  for  the  problem  of  data  association  in  multitarget  tracking.  It  has  also 
been  pointed  out  that  the  DFS  procedure  is  an  efficient  algorithm  for  finding  solutions 
for  the  problem  of  exhaustive  search.  In  this  subsection,  a  specialized  depth-first  search 
DFS  algorithm  for  the  generation  of  data  association  hypotheses  will  be  proposed.  Before 
the  DFS  algorithm  is  given,  it  would  be  worthwhile  to  observe  some  differences  between 
the  problem  of  exhaustive  search  and  that  of  data  association. 

Usually,  there  are  no  solutions  that  are  known  in  advance  in  a  general  exhaustive 
search  problem.  However,  in  the  problem  of  data  association,  a  solution,  which  is  always 
known  in  advance,  is  (0.0 . 0).  The  remaining  solutions  ran  be  generated  system¬ 

atically  from  various  valid  combinations  of  non-zero  values  of  the  elements.  In  order 
to  illustrate  how  this  idea  works,  consider  again  the  example  in  the  previous  sections. 
There  arc  2  targets  and  3  measurements,  and  one  of  the  3  measurements  is  shared  by 
the  2  targets.  A'|,  A'2,  and  X?,  denote  the  3  measurements  and  the  corresponding  Z/s 
are: 

Z\  —  {0, 1};  Z,  =  {0,1,2};  Z3  =  {0,2}. 

The  starting  solution  is  (0,0,0).  The  next  solution  we  look  at  is  (1,0,0),  obtained  by 
setting  A']  -  1.  Since  1  is  already  in  the  previous  solution,  the  only  non-zero  value  A'2 
may  have  is  2.  Therefore,  after  (1,0,0),  the  solution  (1,2,0)  is  obtained.  Now,  all  non-zero 
values  in  Zi  which  can  been  used  to  generate  new  solutions  have  been  used  up.  Reset 
A'2  to  0  and  go  on  to  check  if  there  is  any  non-zero  value  in  Z-{  which  will  lead  us  to  a 
new  solution.  In  this  case,  by  setting  A3  to  2.  another  solution,  (1,0.2),  is  found.  After 
all  non-zero  values  in  both  Zi  and  Z3  are  used  up.  A'i  is  reset  to  0  and  is  kept  as  0  for 
the  remainder  of  the  search  process.  At  this  stage,  A'2  1  will  not  cause  any  conflict. 

Hence,  a  new  solution,  (0,1,0)  is  generated.  This  process  can  go  on  until  all  solutions  arc 
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found.  In  this  particular  example,  the  complete  set  of  solutions  is 


£o  (0,0.0).  <f,  (1.0.0).  £,  (1,2,0). 
£■>,  --  (1,0.2),  S4  -  (0.1.0).  £,  -  (0,1.2), 
U  --  (0,2.0),  £■,  --  (0.0,2). 


For  a  better  view  of  the  solution  generation  process,  all  solutions  can  be  arranged 
into  a  tree,  which  is  called  hypotheses  tree,  as  shown  in  Fig.  2.  Each  node  in  the  tree 


Level  0 


Level  1 


Level  2 


Figure  2  The  graph  arrangement,  of  £i 

stands  for  a  solution  or  a  data  association  hypothesis.  The  root  of  the  tree  is  £0,  which 
implies  that  all  3  measurements  are  hypothesized  to  be  associated  with  clutter.  Each 
node  at  level  /  is  associated  with  i  non-zero  variables  in  the  corresponding  3-tuplc.  This 
also  implies  that  i  out  of  3  measurements  are  hypothesized  to  be  associated  with  i  out 
of  2  targets.  In  Fig.  2.  the  height  of  the  tree  is  2.  In  general,  the  height  of  the  tree  is 
min(n.m)  since  i  must  be  less  than  both  n  and  m. 

As  pointed  out  above,  the  number  of  non-zero  variables  in  a  3-tuple  indicates  the 
level  at  which  the  corresponding  node  is  located.  However,  the  locations  of  the  nodes 
at  the  same  level  are  decided  by  which  variables  are  non-zero  and  their  values.  Suppose 
that  there  are  two  3-tuples  £r  and  £y  and  that  both  of  the  two  3-tuples  have  2  non-zero 
variables.  In  £T,  variables  XJt  and  Xn  are  non-zero  and  in  £y ,  variables  and  Y).  are 
|  non-zero.  £z  is  to  the  left  side  of  £y  if  either 

\ 
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•  A',,  Yj .  and  Xj2  ) when  j{  /,*  and  /.»  jY  or 

•  j i  j;  and  j2  -  ./,*  ■ 

For  example.  £■>  (=  (1,2,0))  and  £>  (-  (1.0.2))  have  2  non-zero  variables.  By  letting  £2 
be  £r  and  £3  be  £y,  we  have  A',  1.  X2  2,  Vi  1.  and  Y>  -  2.  Since  j\  j,*  1  and 

2  "  ./•„>  '  j2*  —  3,  £2  is  at  the  left  side  of  £;{,  as  shown  in  Fig.  2.  Therefore,  the  information 
about  a  data  association  hypothesis  is  contained  in  the  non-zero  variables  in  the  3-tuple. 
This  is  also  true  in  the  case  of  n  targets  and  m  measurements. 

Because  of  the  above  distinguish  features  of  the  data  association  hypotheses,  only 
non-zero  variables  in  a  hypothesis  are  stored  in  a  stack  of  integer  pairs  in  the  imple¬ 
mentation  of  the  DFS  algorithm.  For  example.  £4  (-  (0.1,0))  is  represented  by  ((2,1)). 
The  first  integer  in  the  pair  indicates  which  variable  is  non-zero  in  £4  and  the  second 
integer  is  the  value  of  that  variable.  In  general,  any  m-tuple  with  L  11011-zero  variables 
is  represented  by 

0 '.,*„) 

(j2 1  *£jt) 

where  Xu ,  X}2,  . . .,  and  X]L  are  non-zero  and  L  is  the  pointer  of  the  stack  X  which  is 
pointing  at  the  last  pair  of  integers  in  the  stack.  The  flowchart  of  the  DFS  algorithm  for 
the  generation  of  the  data  association  hypotheses  is  shown  in  Fig.  3.  The  proposed  DFS 
algorithm  is  designated  to  visit  every  node  in  a  hypothesis  tree.  Similar  as  in  Fig.  1, 
GET  X EXT (Zj,  Xj)  in  Fig.  3  is  a  logical  function.  It  is  true  only  if  X}  t  Z;  and  Ar;  is 
not  equal  to  A”|,  A'2,  . . .,  and  XjL .  Otherwise,  it  is  false. 
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Figure  3  The  flowchart  of  DFS  procedure  for  data  association 
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3.3  Computation  of  the  Conditional  Probabilities 


After  each  data  association  hypothesis  is  generated,  the  conditional  probability  of  the 
hypothesis  is  computed  by  using  the  formula  in  (2).  However,  if  the  data  association  hy¬ 
potheses  are  arranged  in  such  a  order  that,  a  factor  of  a  conditional  probability  is  included 
in  another  conditional  probability,  the  computation  of  the  two  conditional  probabilities 
can  be  saved  by  sharing  the  common  factor.  The  data  association  hypotheses  generated 
by  the  proposed  DFS  algorithm  are  in  such  a  special  order  that  the  common  factors  in 
the  conditional  probabilities  can  be  shared  in  the  most  efficient  way.  If  a  little  more 
complex  example  than  the  one  considered  in  the  previous  sections,  the  advantage  of  the 
common  factors  in  the  computation  of  the  conditional  probabilities  may  be  easily  demon¬ 
strated.  Suppose,  there  are  3  targets  and  4  measurements,  measurement  2  is  inside  the 
intersection  of  the  validation  gates  of  targets  1  and  2.  and  measurement  3  is  inside  the 
intersection  of  the  validation  gates  of  targets  2  and  3.  Conventionally,  this  situation  may 
be  represented  by  a  validation  matrix  fl: 


t 


0 

1 

2 

3 

/ 

1 

1 

0 

°] 

1 

1 

1 

1 

0 

2 

1 

0 

1 

1 

o 

<i 

1 

0 

0 

1  \ 

4 

How'cver.  it  may  also  be  represented  by  the  notations  introduced  in  the  mathematical 
model,  exhaustive  search  with  constraints,  proposed  for  data  association.  In  the  second 
representation,  four  variables,  A"),  A^,  X3,  and  A<,  arc  used  to  denote  the  4  measurements 
and  the  corresponding  Z/s  are: 

Z,  -{(),!};  Z2  =  {0, 1,2};  Z3  -  {0,2.3};  Z4  =  {0,3}. 
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When  the  DFS  algorithm  is  applied  to  this  example,  the  following  sequence  of  solutions 
or  data  association  hypotheses  are  generated 


£0  10.0.0.0).  £,  (1.0. 0,0),  El  (1.2. 0.0), 

£,  (1.2. 3.0),  £,  (1. 2.0.3),  £5-  (1. 0.2.0). 

U  (1.0.2, 3).  E-  (1. 0,3,0),  U  (1,0,0, 3). 

The  partial  hypothesis  tree  for  the  solutions  listed  above  is  shown  in  Fig.  4.  Now,  the 

conditional  probabilities  of  the  data  association  hypotheses  Ei  (l  —  0,1,..., 8)  can  be 
computed  using  (2).  Recalling  that  Z  denotes  the  set  of  all  measurements  received  up 
to  the  current  time  index  as  defined  on  page  7. 

P(Eo\Z)  ---  (Po)\  P(£,|Z)  =  (P0)2P/,  P[£i\Z)  =  PoPiPl 

P{E?,\Z)  ---  PI  PE  Pi  P(E4\Z)  ^  PI  P22  Pp  P(E5iZ)  =  PoPjPj, 

P(£g  Z)  P/P/P3,  P(E7\Z)  ---  PoPlPl  P(Es\Z)  =  PoPlPl 

It  is  not  difficult  to  identify  that  P/  is  the  common  factor  \ised  in  the  computation  of 
P(£i\Z)  (/  --  1,2...., 8).  However,  the  common  factors,  P/P|  and  P/P2,  need  a  closer 
examination.  P/ P?  is  used  in  the  computation  of  P(£t\Z)  (l  =  2,3,4),  while  P/P32  is 
used  in  the  computation  of  P(£s|Z)  and  P(£6|Z). 

If  P/’s,  P0,  (Po)2,  and  (Po)3  are  computed  in  advance  for  this  example  and  the  three 
common  factors  mentioned  above  are  used  in  the  computation  of  P(Ei\Z)  (I  0, 1. ...  ,8), 
the  number  of  multiplications  required  in  the  computation  of  P(Ei\Z)  s  may  be  obtained. 
For  example,  there  is  no  multiplication  required  to  obtain  P(Eo\Z).  To  calculate  P(E\  |Z), 
the  product  of  (Po)2  and  P/  is  required.  In  the  computation  of  P(£2\Z),  the  product  of 
P0,  P/ ,  and  P22  is  required.  Therefore,  P{£  1  jZ)  and  P{E2\Z)  are  computable  with  one  and 
two  multiplications.  Since  the  product,  P/  P22,  is  obtained  in  the  computation  of  P(£2\Z), 
there  is  only  one  multiplication  required  to  compute  of  P(<f3|Z)  (  =  P/P2  P33).  The  number 
of  multiplications  required  in  the  computation  of  the  conditional  probabilities.  P(£{|Z)’s. 
of  the  data  association  hypotheses  at  different  levels  in  the  hypothesis  tree  is  listed  in 
Table- 1. 
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In  general,  the  computation  for  Pj .  Pa,  (Pa)2 . and  is  called  pre¬ 

processing,  which  facilitates  the  use  of  common  factors  in  the  computation  of  the  condi¬ 
tional  probabilities  of  the  data  association  hypotheses.  The  details  of  the  multiplications 
required  in  the  computation  of  the  conditional  probabilities  of  the  data  association  hy¬ 
potheses  at  each  level  of  the  hypothesis  tree  is  summarized  in  Table-2. 


Let  N D,  denote  the  number  of  nodes  at  level  i  and  M  be  the  total  number  of 
multiplications  required  in  the  computation  of  P(£/jZ)’s.  From  Table-2,  it  can  be  inferred 
that 

min(n.m)  -  1 

XI  ND\  ,  2  £  ND,  - NDminnm 

min  ( n .  ra ) 

<  2  £  (ii) 

i-  o 

(11)  shows  that  M  is  less  than  twice  the  number  of  data  association  hypotheses  if  pre¬ 
processing  and  common  factors  are  used.  This  enables  11s  to  analyze  the  computational 
complexity  of  the  DFS  algorithm  in  terms  of  the  number  of  data  association  hypotheses 
or  the  number  of  nodes  in  the  hypothesis  tree. 


Level  0 


Level  1 


Level  2 


Level  3 


Figure  4  The  partial  hypothesis  tree  for  the  example  given  in  this  section 
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Table  1  Number  of  multiplications  required  in  the  computat.i  .11  of  P(£i  Z) 's 


(Conditional  Probabilities  Number  of  Multiplications  i 

_P(£»:Z)^  _  0 

P(£i \Z)  _  _  1 

P{£2  Z)^P{£,\Z),  P{£7\Z),  P{£h\Z)  2  __ 

P{£vZ),  P(£4'Z).  P(£b\Z)  1 


Table  2  Number  of  multiplications  required  in  the  computation  of  P(£i\Z)' s  in  genera 


Hypotheses  at  Level  i  Number  of  Multiplications 


min(n,m)  1  j  2 

min(u,m)  j  1 
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3.4  Computational  Complexity 


With  the  preprocessing  discussed  in  the  previous  subsection,  the  computational  com¬ 
plexity  analysis  is  simplified  as  counting  the  number  of  data  association  hypotheses  or 
the  number  nodes  in  the  hypothesis  tree.  In  general,  as  we  can  expect,  it  is  impossible 
to  count  the  number  of  the  nodes  in  the  hypothesis  tree  for  any  situation.  In  this  sub¬ 
section.  the  computational  complexity  is  analyzed  in  two  cases,  the  worst  case  and  the 
average  case.  Suppose  that  there  are  n  targets  and  rn  measurements.  According  to  the 
mathematical  model  for  data  association.  rn  variables.  X}  (j  =  1,2,.  .  .  . m),  are  used  to 
denote  the  rn  measurements.  Each  variable  X}  is  defined  on  a  set  Z;.  In  the  worst  case. 
Zj  assumes  the  following  form: 

Zj  ■---  {0. 1, 2, ....  n}  for  j  -  1,2  (12) 

In  (12),  each  Zj  contains  n  non-zero  values.  In  the  average  case,  however,  the  number 
of  non-zero  values  in  each  Z3  is  much  less  than  the  total  number,  n,  of  targets.  In  the 
following,  the  worst  case  analysis  is  considered  first.  Then,  the  average  case  analysis  is 
discussed. 

3.4.1  The  Worst  Case  Analysis 

Let  S(n.rrt)  denote  the  number  of  data  association  hypotheses  or  the  number  of  nodes 
in  the  hypothesis  tree  in  the  worst  case.  Suppose,  n  =  1.  In  this  case,  either  this  target 
is  hypothesized  to  be  associated  with  one  of  the  m  measurements  or  all  m  measurements 
are  hypothesized  to  originate  from  clutter.  Therefore,  5(1,  m)  =  rn  +  1.  Similarly,  when 
rn  1,  the  total  number  of  data  association  hypotheses,  5(n,  1),  is  equal  to  n  +  1.  In 
general,  the  DFS  algorithm  starts  at  the  root  of  a  hypothesis  tree,  where  all  variables  are 
zero.  Then.  A"i  is  set  to  1.  The  total  number  of  valid  combinations  of  non-zero  values  of 
the  variables  X}  ( j  ~  2, 3, . . . .  in)  is  S(n  -  1,  rn  -  1).  Since  Xj  —  1,  the  remaining  m  1 
variables  X}  (j  —  2.3, . . .  ,  rn)  can  only  have  n  1  non-zero  values,  i.e.,  2,  3. ...  ,  n.  After 
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all  combinations  of  non-zero  values  of  the  variables  Xj  (j  -  2.3 . rn)  are  generated, 

A’ i  is  set  to  2.  As  just  discussed,  the  total  number  of  valid  combinations  of  the  variables 
Xj  (j  -  2,3,..  .  .  rn)  is  S(n  l.m  1).  This  process  can  go  on  until  all  non-zero  values 
of  A'i  are  used.  Since  there  are  n  different  non-zero  values  for  variable  X\,  the  number 
of  data  association  hypotheses  in  which  Aj  is  non-zero  is  rtS(n  -  l.m  —  l).  Once  all  the 
data  association  hypotheses  with  X]  /  0  are  generated,  the  DFS  algorithm  goes  back 
to  the  root  of  the  hypothesis  tree  and  resets  A't  to  0.  The  remaining  data  association 
hypotheses  are  then  generated  based  on  rn  -  1  variables  A",  (j  =  2,3,...,  m)  and  n  non¬ 
zero  values  for  each  variable.  Therefore,  the  number  of  data  association  hypotheses  with 
A’i  -  0  is  5(n,  rn  —  1).  It  is  obvious  now  that  the  number  of  data  association  hypotheses 
or  the  number  nodes  in  a  hypothesis  tree  in  the  worst  case  can  be  represented  by  a 
difference  equation, 

S(n.m)  -  nS(n-l,m  1)  4  S(n.m  -  l).  (13) 

with  boundary  conditions. 

5(1,  m)  -  m  +  1, 

S(n,  1)  -  n  4  1. 


Since  (13)  is  not  a  constant  coefficient  difference  equation,  its  solution  is  not  very 
easy  to  get.  To  avoid  the  trouble,  return  to  the  structure  of  the  hypothesis  tree.  Recalling 
that  the  nodes  at  level  i  represent  the  hypotheses  that  i  out  of  m  measurements  are  from 
i  out  of  n  targets.  There  are  C'm  choices  to  get  any  i  out.  of  rn  measurements,  where 


C 


m 


ml 


for  1  <  i  <  rn 


1  for  i  0  . 


V 


Among  n  targets,  there  are  A'n  different  ways  of  associating  the  i  measurements 
with  the  i  out  n  targets,  where 


[  n(n  1)  •  •  •  (n  -  i  4  1)  for  1  <  i  <  n 

|  1  for  i  —  0  . 
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So.  the  number  of  nodes  in  the  tree  is 


min(n,m) 

S(n.m)  Y.  A‘>  Cm-  (14) 

i  —0 

It  can  be  shown  that,  the  solution  satisfies  (13).  From  the  solution,  it  is  apparent  that 
the  computational  cost  is  highly  increasing  with  the  increase  of  n  and  m.  In  Table-3, 
some  S{n,rn)'s  are  listed  with  different  numbers  of  targets  and  measurements.  From  the 
listed  values,  we  can  see  that  even  when  the  number  of  targets  is  as  small  as  six,  the 
number  of  nodes  in  the  hypothesis  tree  could  be  greater  than  one  million  in  the  worst 
case. 


Table  3  Sample  values  of  S(n,m) 


n  i  m 

I  _ 

— 

S(n,m) 

n 

rn 

S(n,m) 

j 

2  j  4 

21 

4 

* 

3,393  : 

3  6 

u _ 

229 

_ 

6 

12 

1,442,173  |i 

In  practical  situations,  however,  when  the  densities  of  target  and  clutter  are  not  very 
high,  the  number  of  non-zero  values  in  each  Z2  is  much  less  than  that  of  measurements, 
m.  Consequently,  the  number  of  hypotheses  is  much  less  than  S(n,  m).  One  can  derive  a 
tighter  upper  bound  by  introducing  additional  assumptions  to  modify  (13).  For  example, 
the  number  of  non-values  in  each  Z}  is  less  than  or  equal  to  two.  which  means  that  any 
measurement  cannot  be  shared  by  more  than  two  targets.  With  this  assumption,  (13) 
may  be  changed  to 


5(n,m)  =  2 S(n  —  l,m  -  1)  f  S(n,m  -  l).  (15) 

In  next  subsection,  a  simple  practical  upper  bound  for  the  number  of  data  association 
hypotheses  is  derived  in  the  average  case  where  the  number  of  non-zero  values  for  each 
variable  A';  is  much  less  than  that  of  measurements,  m. 
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3.4.2  The  Average  Case  Analysis 


Commonly,  it  is  assumed  that  there  are  on  the  average  two  measurements  inside  the 
validation  gate  of  each  target.  In  other  words,  the  number  of  non-zero  values  in  each 
Z}  is  less  than  or  equal  to  two  under  this  assumption.  For  each  target,  there  are  three 
different  associations  with  the  two  measurements.  That  is.  the  target  is  hypothesized  to 
be  associated  with  one  of  the  two  measurements  and  the  target  is  not  associated  with  any 
one  of  the  two  measurements.  Therefore,  on  the  average,  the  number  of  data  association 
hypotheses  is  bounded  by 

n 

S*(n,m)  =-  3  *  3  *  •  •  -  *  3  -  3".  (16) 

In  the  example  discussed  in  Section  2.  there  are  2  targets  and  3  measurements  and 
there  are  2  measurements  inside  the  validation  gate  of  each  target.  The  number  of  data 
association  hypotheses  obtained  from  the  worst  case  analysis  is  13  and  the  upper  bound 
obtained  from  ( 16)  is  9.  The  actual  number  of  data  association  hypotheses  in  this  case  is 
8.  The  upper  bound  given  in  (16)  can  be  easily  extended  to  the  case  where  the  number 
of  measurements  inside  the  validation  gate  of  each  target  is  ml .  That  is, 

n 

S*(n,m)  ~  JI(l  r  rn‘).  (17) 

i- 1 
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4  Direct  Computation  of  3' 


As  discussed  in  the  previous  section,  the  computational  cost  for  the  generation  of  the  data 
association  hypotheses  is  very  high  when  ti .  the  number  of  targets  in  a  cluster,  and  m, 
the  number  of  measurements,  are  relatively  large.  To  further  speed  up  the  computation 
of  the  a  posteriori  probability  /i),  an  algorithm  which  applies  to  certain  special  situation 
is  proposed  in  this  section.  This  procedure  does  not  requires  the  generation  of  the 
data  association  hypotheses.  This  algorithm  is  applicable  when  the  density  of  targets 
is  moderate.  Specifically,  the  formulae  for  computing  j3l  directly  are  developed  for  the 
cases  when  n  equals  1,  2.  3,  or  4. 


When  a  cluster  is  composed  of  only  one  target  the  PDAF  is  applied.  Using  the 
notations  in  (3)  and  (4),  we  have 

-  Pj  for  j  1,2, ... ,  m,  (18) 

4  -  Po.  (19) 

Every  3‘  is  normalized  by 

m  m 

c  a-£p;  £4  (20) 

/-i  l=o 

The  normalized  is  the  one  occurring  in  the  PDAF  presented  in  [1|  with  Pc  =  1,  where 


Pc  is  the  probability  for  the  correct  measurement  to  fall  inside  the  validation  gate. 

For  the  example  given  in  Section  2,  there  are  2  closely  spaced  targets  and  3  measure¬ 
ments  are  received.  One  of  the  3  measurements  is  inside  the  intersection  of  the  validation 
gates  of  the  two  targets.  Therefore,  these  two  targets  are  in  the  same  cluster.  In  this 
case,  each  3l}  is  a  summation  of  the  conditional  probabilities  of  some  data  association 
hypotheses.  For  example, 

0\  =  P(£(ih)\Z)  + P{£{tl2)\Z)  f  P(£((h)\z). 

If  the  notations  in  (3)  and  (4)  are  used,  the  above  equation  can  be  simplified  as, 

P/P0  t  P,‘P22  +  Pi  Pi 
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p,1(p<>  ■  Pi  Pi). 


(21) 

In  (21).  there  is  no  reference  to  the  conditional  probabilities  of  data  association  hypothe¬ 
ses  n)5  il2,  and  n>.  it*  is  directly  computed  from  Pj  and  P„.  In  general,  for  2  targets  and 
7,1  measurements,  d)  can  be  obtained  by  the  following  equations.  For  j  —  1,2, 

/,  -  1,2.  for  i  -  1.2.  and  /|  /  /_>, 

ft  P‘/(Ptr  P ‘2)  (22) 

ft  PoP h  (23) 

where 

m 

-  P0  t  £  P;(  for  t  —  1,2.  (24) 

j^i 

In  (22),  the  term  P  1 P^ 2  is  dropped  because  Pj‘ P*2  is  the  conditional  probability  for 
measurement  j  to  originate  from  target,  t\  and  I21  which  violates  the  two  restrictions  [3] 
given  in  Section  2. 


Similarly,  for  3  targets  and  m  measurements,  3)  can  be  computed  by  the  family 
of  equations  given  below.  For  j  -  1. 2. . . . ,  m,  /,  -=  1. 2. 3.  for  1  =  1.2.3,  and  tp  £  tq  if 
P  *  <7, 


m 


ft 

--  P,,,[(P'2 

P,‘2)(P(’  -  P/3)  -  £  P/2P/3!, 

i  1 

/  r  1 

(25) 

ft 

=  Po(PhPh 

‘  7  J 

Eft'*). 

(26) 

When  the  number  of  targets  in  a  cluster  is  increased  to  4,  the  set  of  equations  for 
3]  becomes  much  more  complex.  For  j  =  l,2,...,m,  t,  1,2, 3,4,  for  i  ■=  1.2, 3,4,  and 
tp  r  tq  if  P  3  q, 


25 


(28) 


y'1 

■  'n 


p0( n  r1'  Phl  V  p/1  p/«  p'1  V  p/*p/« 


/  - 1 


i 


P‘<  V  pt‘*p/»  -  2  V  P/2P/SP/4) 


/  i 


/  - 1 


When  the  number  of  targets  is  larger  than  4,  the  equations  for  computing  /P  can 
be  derived  with  careful  consideration.  However,  the  complicated  equations  themselves 
may  require  a  lot  of  memory.  In  a  small  system,  excessive  requirement  for  memory  may 
not  be  suitable.  Therefore,  in  this  section,  the  equations  for  direct  computation  of  /?' 
are  derived  for  the  cases  in  which  the  number  of  targets  in  a  cluster  is  less  than  or  equal 
to  4.  The  normalizing  constant  c  in  all  cases  is  obtained  from 

TO 

c  -  £#  (29) 

/  o 

In  to  order  to  compare  with  the  DFS  algorithm,  the  number  of  multiplications 
required  in  the  computation  of  li‘  in  the  algorithm  proposed  here  is  estimated  from 
the  above  equations  in  all  cases.  Since  the  computational  costs  for  obtaining  /?]  in  the 
two  algorithms  are  the  same  when  there  is  only  one  target  in  a  cluster,  the  following 
discussion  is  focused  on  the  cases  in  which  the  size  of  a  cluster  of  targets  is  larger  than 
1.  In  the  worst  case,  in  the  case  of  2  targets  and  m  measurements,  one  multiplication  is 
required  by  (22)  and  one  by  (23).  Therefore,  the  total  number  of  multiplications  required 
to  compute  all  tfj’s  in  this  case  is  2(m  -  1).  If  the  DFS  algorithm  is  applied,  the  total 
number  of  multiplications  can  be  obtained  by  using  (11).  Suppose,  rn  >  2.  Then. 

A/  -  /Ijt’i  f  A\Cl 

—  m(m  +  1) 

>  2(m  t  1)  (30) 

where  A'2C'm  is  the  number  of  nodes  at  level  i  in  the  hypothesis  tree  as  discussed  in  the 
complexity  analysis  in  the  previous  section.  Similar  comparisons  can  be  made  between 
the  DFS  algorithm  and  the  direct  computation  formulae  for  (3‘  in  the  case  of  3  and  4 
targets,  as  shown  in  Table-4.  In  Table-4,  DC  and  DFS  stand  for  direct  computation 
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and  depth-first  search  algorithms,  respectively.  The  results  in  Table-4  are  based  on  the 
assumption  that  in  is  larger  than  or  equal  to  n.  Generally  speaking,  the  number  of 
multiplications  required  in  the  direct  computation  of  J*  is  less  than  that  required  in  the 
DFS  algorithm  as  shown  in  Table-4.  However,  this  is  not  true  when  n  and  m  are  equal  to 
4.  Therefore,  it  is  quite  possible  that  the  number  of  multiplications  required  in  the  direct 
computation  algorithm  could  be  larger  than  that  required  in  the  DFS  algorithm  when 
n  :•  4  and  rn  is  almost  equal  to  n.  In  the  implementation,  it  would  be  a  good  idea  to 
combine  the  two  algorithm  together.  For  each  cluster  of  targets,  the  direct  computation 
algorithm  is  applied  if  n,  the  size  of  the  cluster,  is  less  than  or  equal  to  4.  Otherwise, 
the  DFS  is  applied. 


Table  4  Number  of  multiplications  required  in  the  two  algorithms 


:  »■ 

DC 

DFS 

i 

2 

2 (m  -4- 

1) 

m(m  + 

1) 

j 

i  O 

1  ° 

3(mJ  -f-  2 rn  s 

2) 

m(m2  +  3m  - 

1) 

!"  4  " 

1 

4(5m2  -f  G rn  ■* 

G)_ 

rn(ni'  -  2 m:  rn  - 

2) 
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5  Approximation  of  .V 


Unlike  the  two  algorithms  which  have  boon  developed  in  the  previous  sections  to  compute 
3‘  efficiently,  in  this  section,  an  algorithm  is  proposed  to  compute  .ij  approximately 
without  computing  the  conditional  probabilities  of  the  data  association  hypotheses  in 
the  JPDAF.  In  the  original  JPDAF,  the  interference  from  all  the  targets  in  the  same 
cluster  is  considered  in  the  computation  of  .  This  results  in  high  computational  cost 
for  the  data  association.  To  reduce  the  computational  cost  for  the  data  association, 
only  the  interference  from  the  neighbors  of  target  /  is  considered  in  the  approximation 
of  ;V ■  A  target  is  a  neighbor  of  target  t  if  there  is  at  least  one  measurement  inside 
the  intersection  of  the  validation  gates  of  the  two  targets.  It  is  not  necessary  that  two 
targets  in  the  same  cluster  are  neighbors.  For  example,  targets  1,  2.  and  3  are  in  the  same 
cluster  if  targets  1  and  2  are  neighbors  of  target  3.  However,  if  there  is  no  measurement 
inside  the  intersection  of  the  validation  gates  of  targets  1  and  2,  targets  1  and  2  are  not 
neighbors. 

In  the  direct  computation  algorithm,  each  (i)  is  computed  as  following. 

3]  -  P‘ Fj  for  t  -  l,2,...,n,  and  j  =  0,1,2 . rn.  (31) 

In  (31).  Fj  is  the  interference  from  the  other  targets  in  the  same  cluster.  Comparing 
(31)  with  (18),  the  interference  Fj  may  be  considered  as  a  weight  on  Pj.  To  see  what 
effect  of  the  weight  Fj  has  on  the  value  of  /?),  let  us  consider  the  example  in  Section  2. 
In  that  example,  there  are  2  targets  and  3  measurements.  Among  the  3  measurements, 
measurement  2  is  inside  the  intersection  of  the  validation  gates  of  the  two  targets.  Using 
(22)  and  (23).  we  have 


Po  ( Po  *  Pi  i 

Pi)  -  FoFj, 

(32) 

P}(Po  +  P2  + 

Pi)  =  F/F/, 

(33) 

P2'(P0  +  Pi) 

P2‘ f2‘. 

(34) 

From  the  above  equations,  it  is  obvious  that  Fj  is  less  than  both  F0'  and  F,1.  This  means 


that  more  weight  is  put  on  measurement  1  than  on  measurement  2  since  measurement 
2  is  inside  the  intersection  of  the  validation  gates  of  the  two  targets.  In  other  words, 
a  measurement  is  less  likely  to  come  from  target  t  if  the  measurement  is  inside  the 
intersection  of  the  validation  gates  of  target  t  and  other  targets.  If  Fj's  arc  equal  for 

j  0, 1 . m,  then  there  is  either  no  interference  from  any  other  target  or  the  effect 

of  the  interferences  on  the  values  of  /?j’s  are  equal.  This  happens  when  there  is  no 
intersection  between  the  validation  gates  of  target  t  and  the  other  targets.  The  value  of 
,3*  obtained  from  the  joint  data  association  with  the  other  targets  is  the  same  as  that  of 
3l-  from  the  data  association  for  target  t  alone.  Generally,  if  targets  are  not  grouped  into 
clusters,  the  value  of  3l  for  each  target  and  each  measurement  is  the  same  as  that  of  3‘ 
after  targets  are  grouped  into  clusters.  However,  much  more  computation  is  needed  to 
compute  3j  if  targets  are  not  grouped  into  clusters.  So,  it  is  very  important  that  targets 
are  grouped  into  clusters  before  considering  data  association. 

Even  after  clustering,  the  computation  to  obtain  F‘  could  be  still  very  costly  if  n, 
the  number  of  targets  in  a  cluster,  and  m,  the  number  of  measurements,  are  relatively 
large  as  discussed  in  the  previous  sections.  To  reduce  the  computational  cost,  in  the 
approximation  of  3),  only  the  interference  from  the  neighbors  of  target  t  is  considered 
in  the  computation  of  Fj .  Let  Dl  denote  the  set  of  neighbors  of  target  t.  That  is,  a 
target  belongs  to  J3(  if  there  is  at  least  one  measurement  inside  the  intersection  of  the 
validation  gates  of  this  target  and  target  t.  Then,  the  computation  of  F‘  is  simplified  as 
shown  below. 

m 

F‘  po  +  E  E  Pi  for  j-  1,2, - m.  (35) 

i€0'  l  =  I 
l  *  J 
m 

EE  PI-  (36) 

i</n  i=i 

Substitute  Fj  into  (31),  to  obtain 

-  PjFj1  for  t  --  1,2 ,...,n,  and  j  -0,1,2, _ m.  (37) 

Finally.  3]  is  normalized  by  the  constant  c  in  (29). 
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From  the  manner  in  which  F‘  is  computed,  this  algorithm  is  actually  a  approx¬ 
imation  of  the  direct  computation  algorithm  presented  in  the  previous  section.  When 
n  -  2.  (37)  is  equivalent  to  (22)  and  (23).  Comparing  with  PDAF  ! 2 , ,  the  computational 
cost,  of  rhis  approximation  algorithm  is  increased  slightly  because  of  computing  F!.  Fur¬ 
thermore.  it  is  expected  that  the  performance  of  the  approximation  algorithm  could  be 
close  to  that  of  the  JPDAF  because  the  interference  from  neighboring  targets  of  target  / 
is  considered  in  the  computation  of  0l.  The  computer  simulation  of  the  approximation 
algorithm  is  presented  in  the  next  section. 
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6  Computer  Simulation 


In  the  simulation,  the  dynamic  models  for  the  targets  have  been  digitized  using  the 
sampling  period  T  normalized  to  l.s  and  the  state  vectors  have  been  represented  in 
2-dimensional  Cartesian  coordinates  (or  in  A'  1"  plane).  Furthermore,  only  position 
measurements  have  been  assumed  to  be  available.  The  surveillance  region  used  in  the 
simulation  is  a  35/cm  by  35 krn  square  and  the  initial  positions  and  velocities  of  16  targets 
in  2-dimensional  Cartesian  coordinates  are  given  in  Table-5.  The  performance  of  the 
three  algorithms  proposed  are  tested  in  two  separate  cases  in  the  simulation.  In  the  first 
case,  suppose  that  the  dynamic  characteristics  of  all  targets  are  known,  i.e.,  every  target 
moves  at  a  constant  velocity.  In  the  second  case,  however,  some  targets  could  make  a 
maneuver  at  any  time.  In  the  simulation,  target  10  and  11  are  capable  of  making  a  turn 
with  a  centripetal  acceleration.  The  maneuver  parameters  of  the  two  targets  are  given 
in  Table-6. 

Since  every  target  is  nonmaneuvering  in  the  first  case,  the  generic  target  dynamic 
model  has  the  following  form: 


x(fc  +1)  =  Fx{k)  •+  Gw(/c) 

(38) 

z  (k)  -  Hx(k)  +  v(/c), 

(39) 

where 

x(^)  (  x{k)  x(k)  y(k)  y{k)  )  •  w (k)  -  (  i v,(k)  w2(k)  'j  ■ 

E{w(k.) }  -  0,  £{w(fc)w'(j)}  -  Q6(k  j), 
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l  0 

r  J 
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1  0  0  0 


0  0  10 


E{v(k)}  -  0.  E{v{k)v'(j)}  R6(k-j). 

For  maneuvering  targets,  a  higher  order  of  dynamic  model  is  required  to  describe 
the  characteristics  of  these  targets.  In  the  simulation,  the  modified  Singer  s  model  de¬ 
veloped  in  1 1 2 j  is  used  for  all  targets  in  the  second  case.  The  state  and  observation 

equations  in  the  modified  Singer’s  model  are  described  by 

xm(k  -f  1)  -  FmXm(k)  r  u(k)  (40) 

z(k)  -  HmXm(k)  \-x(k).  (41) 

In  (40)  and  (41),  xm(fc),  Fm,  and  Hm  are  defined  as 

xm(Ar)  (  x(k)  x(k)  x(k)  y(k)  y(k )  y(k)  )  , 


f(T)  0 
0  f(T) 


1  0  0  0  0  0 


0  0  0  1  0  0 


1  T  T2i  2 
0  1  T 


0  0  1 


The  statistical  properties  of  v(k)  in  (41)  are  identical  to  those  of  v(k)  in  (39).  In  (40)  and 
(41),  the  two  noise  terms,  u [k)  and  v(/c)  are  still  assumed  to  be  mutually  independent. 
The  state  noise  u(fc)  is  a  zero-mean  Gaussian  white  random  process  with  covariance 
matrix  Qm{k.)  of  the  following  form: 


Qm(k) 


ox(k)q(T)  0 
0  oy(k)q(T) 


r5  -  20  t*/  8  r3/  e 

TV  8  T3/ 3  T2)  2  , 

r3/ 6  r-,/2  t 
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where'  <>  is  the  maneuver  correlation  coefficient,  and  of(k)  and  oy(/c)  are.  respectively, 
the  .V  and  1'  components  of  the  maneuver  magnitude  variance.  Both  a T(k)  and  oy(k) 
are  adaptively  estimated  in  the  simulation. 

The  correct  returns  from  a  target  are  generated  by  adding  noise  to  the  computed 
true  position  of  the  target.  The  standard  deviation  of  the  measurement  noise  has  been 
selected  as  y/Ra  --  0.15 km  for  both  the  X  and  Y  vomponents.  The  correct  return  would 
pass  a  detector  with  probability  of  detection  Pp  0.99.  The  clutter  is  generated  uni¬ 
formly  over  the  whole  surveillance  region.  The  total  number  of  clutter  returns  observed 
in  the  region  is  a  Poisson  random  number.  The  density  of  the  clutter.  A,  is  selected  to 
be  0.05  km2.  This  has  given  an  average  of  0.5-2  clutter  per  gate  in  the  examples  given 
below.  The  threshold  g 2  for  the  validation  gate  is  set  to  17.0. 

The  efficiency  of  the  three  algorithms  is  estimated  in  terms  of  CPU  time  used  in 
the  data  association  in  a  single  sample  run.  The  simulation  is  performed  on  a  VAX  8550. 
The  exact  CPU  time  for  each  example  could  vary  from  time  to  time,  depending  upon 
the  load  of  the  computer.  The  robustness  of  the  three  algorithm  is  evaluated  using  the 
success  rate  in  tracking  both  maneuvering  and  nonmaneuvering  targets  in  clutter.  The 
su<  css  rate  means  the  number  of  sample  runs  finished  without  lost  tracking  of  a  single 
target  in  100  different  sample  runs. 

In  the  simulation,  three  examples  are  used  to  track  5.  9,  and  16  targets  with  different 
dynamic  models  given  in  (38)  and  (40).  The  initial  parameters  of  the  5  targets  in  the 
first  example  arc  taken  from  targets  4,  6,  7,  11.  and  13  and  the  initial  parameters  of  the 
„  targets  in  the  second  example  are  taken  from  targets  1,  2,  4,  6,  7,  8,  10.  11.  and  16  in 
Table-5.  In  the  third  example,  all  targets  listed  in  Table-5  are  used.  When  all  targets 
are  supposed  to  be  nonmaneuvering,  the  three  examples  are  shown  in  Figs.  5.  7,  and  9. 
If  the  modified  Singer's  mod'd  in  (40)  is  used  for  all  targets,  the  trajectories  of  all  targets 
in  the  three  examples  are  shown  in  Figs.  6,  8,  and  10. 

The  CPU  times  used  in  the  data  association  in  one  sample  run  in  the  three  algo- 


rit bins  arc'  listed  in  Table-7.  In  Table-7.  DFS.  AC.  and  DC  stand  for  depth-first  search, 
approximate  computation,  and  direct  computation  algorithms,  respectively.  In  the  first, 
example,  as  shown  in  Figs.  5  and  0.  the  largest  size  of  cluster  of  targets  is  3.  But  most 
of  the  time  during  simulation,  the  sizes  of  clusters  are  1  and.  occasionally.  2.  When 
there  is  one  target  in  a  cluster,  a  PDAF  is  applied  for  tracking  this  target.  In  the  three 
algorithms,  the  equations  used  to  compute  /Jj  are  the  same  when  the  size  of  a  cluster  is  1. 
As  a  result,  in  this  example,  the  CPU  times  used  in  the  data  association  are  almost  the 
same  in  the  three  algorithms  for  tracking  nonmaneuvering  and  maneuvering  targets.  In 
the  second  example,  as  shown  in  Figs.  7  and  8.  the  total  number  of  targets  is  increased 
to  9  and  the  largest  size  of  cluster  of  targets  is  increased  to  4.  Most  of  the  time  in  the 
simulation,  the  size  of  the  largest  cluster  at  each  time  instant  is  about  2.  Therefore, 
the  increases  of  the  CPU  times  in  the  three  algorithms  are  roughly  proportional  to  the 
increase  of  the  number  of  targets.  However,  in  the  third  example  as  shown  in  Figs.  9  and 
10.  the  size  of  the  largest  cluster  is  frequently  ranging  from  6  to  8,  and  it  is  sometimes  as 
large  as  9.  This  leads  to  a  significant  increase  of  the  CPU  time  used  for  data  association 
in  the  DFS  algorithm.  For  the  approximate  computation  (AC)  algorithm,  the  increase 
of  the  CPU  time  is  proportional  to  the  increase  of  the  number  of  targets  because  of  the 
way  in  which  the  data  association  is  considered  in  this  algorithm.  Although  the  AC 
algorithm  is  much  more  efficient  than  the  DFS  algorithm  when  the  targets  are  nonma¬ 
neuvering,  the  AC  algorithm  is  not  suitable  for  use  to  track  16  targets  when  there  are 
possible  maneuvers.  Since  the  size  of  the  cluster  of  targets  is  frequently  larger  than  4  in 
this  example,  the  direct  computation  (DC)  algorithm  is  not  applicable. 

The  success  rates  of  the  three  algorithms  are  listed  in  Table-8.  When  targets  are 
nonmaneuvering,  target  11,  as  shown  in  Fig.  9,  causes  most  of  the  mistracking.  However, 
when  the  dynamic  model  given  in  (40)  is  used,  the  two  maneuvering  targets  are  the 
major  cause  of  the  mistracking.  Of  course,  the  different  layout  of  the  targets  could  lead 
to  different  success  rates.  Especially,  the  success  rates  might  be  improved  if  either  target 
11  is  removed  or  it  is  moved  to  a  different  position.  However,  the  current  layout  of  the 
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targets  could  help  us  understand  the  limitation  of  the  joint  data  probabilistic  association 
technique.  In  all  examples,  the  available  success  rates  of  the  DFS  and  the  DO  algorithms 
are  the  same,  as  expected,  when  the  density  of  targets  is  not  high.  But  the  DO  algorithm 
is  computationally  more  efficient  than  the  DFS  algorithm.  The  high  success  rate  of  the 
DFS  algorithm  is  accompanied  by  high  computational  cost  in  comparison  with  the  AC 
algorithm.  The  success  rate  of  the  AO  algorithm  in  the  case  of  tracking  16  maneuvering 
targets  is  not  available  in  Table-8  due  to  the  same  reason  given  above.  The  performance  of 
the  DFS  and  the  AC  algorithms  deteriorates  with  the  increase  in  the  number  of  targets. 
This  difficulty  is  caused  by  the  increase  of  the  number  of  measurements  inside  each 
gate.  Therefore,  additional  information  is  needed  to  reduce  the  uncertainty  in  the  data 
association.  Practically,  the  choice  of  the  three  algorithm  in  the  design  of  a  tracking 
system  will  be  dependent  upon  the  computational  capability  and  the  requirement  of 
tracking  accuracy. 
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Table  5  The  initial  positions  and  velocities  of  1G  targets 
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Table  6  The  maneuver  parameters  of  target-10  and  target-11 
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Table  7  CPU  Time  Used  in  Data  Association  in  Seconds  (one  sample  run) 
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Table  8  The  success  rates  of  the  three  algorithms  in  100  runs 
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Figure  5  Tracking  5  nonmaneuvering  targets. 
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Figure  6  Tracking  5  maneuvering  targets. 
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Figure  7  Tracking  9  nonmaneuvering  targets. 


40 


k 


40 


X  (km) 


X 

Clutter 

Estimated  Track 

& 

Target  Return 

- True  Track 

o 

Predicted  Position 

Figure  8  Tracking  9  maneuvering  targets. 
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Figure  9  Tracking  16  nonmaneuvering  targets. 
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Figure  10  Tracking  16  maneuvering  targets. 


7  Conclusion 


In  this  paper,  three  algorithms  have  been  proposed  for  fast  computation  of  the  a  posteri¬ 
ori  probability  in  the  JPDAF.  As  shown  in  the  simulation,  the  CPU  time  is  increased 
drastically  when  the  largest  size  of  the  cluster  of  targets  is  relatively  large.  Therefore,  the 
DFS  algorithm  is  suitable  for  implementation  in  a  tracking  system  and  in  a  ground-based 
surveillance  system  with  a  large  centralized  computational  capability.  Furthermore,  the 
JPDAF  implemented  with  the  DFS  algorithm  could  serve  as  a  standard  for  any  future 
attempt  to  approximate  the  .JPDAF.  The  advantages  of  the  AC  algorithm  proposed  in 
this  paper  are  that  the  computation  of  the  a  posteriori  probability  is  very  efficient 
and  that  the  success  rate  is  relatively  high  in  comparison  with  the  JPDAF  without  any 
approximation.  Especially,  the  AC  algorithm  is  suitable  for  implementation  in  a  multi¬ 
processor  system.  The  DC  algorithm  is  more  efficient  than  the  DFS  algorithm  and  can 
be  implemented  in  a  multiprocessor  system.  However,  right  now,  it  is  only  applicable 
when  the  density  of  targets  is  not  high. 

Although  the  discussion  here  is  focused  on  the  target-oriented  approach,  i.e.,  JPDAF 
3!.  it  is  not  difficult  to  extend  the  DFS  algorithm  to  the  measurement-oriented  approach 
4j.  Furthermore,  the  DFS-based  technique  is  applicable  to  three  dimensional  validation 
matrices  which  occur  in  Markov  models  of  system  parameters  for  maneuvering  targets 
in  clutter  }  10 j  and  in  multiscan  correlation  [4 j . 
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